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Summary 

The efficiency and accuracy of several algorithms 
recently developed for the efficient numerical integration 
of stiff ordinary differential equations are compared. 
The methods examined include two general-purpose 
codes, EPISODE and LSODE, and three codes 
(CHEMEQ, CREKID, and GCKP84) developed specifi- 
cally to integrate chemical kinetic rate equations. The 
codes are applied to two test problems drawn from 
combustion kinetics. The comparisons show that LSODE 
is the fastest code currently available for the integration 
of combustion kinetic rate equations. 

An important finding is that an iterative solution of the 
algebraic energy conservation equation to compute the 
temperature does not result in significant errors. In addi- 
tion, this method can be more efficient than evaluating 
the temperature by integrating its time derivative. 

SigniEcant reductions in computational work are 
realized by updating the rate constants (k=AT^ 
exp( -£■//? 7)) only when the temperature change exceeds 
an amount AT that is problem dependent. An 
approximate expression for the automatic evaluation of 
AT is derived and is shown to result in increased 
efficiency. 

Introduction 

The major problem associated with the simultaneous 
numerical integration of large sets of chemical kinetic 
rate equations is that of stiffness. Although stiffness does 
not have a simple definition (see, e.g., Shampine, refs. 1 
and 2), it is characterized by widely varying time 
constants. For example, in hydrogen-air combustion the 
induction time is of the order of microseconds whereas 
the nitric oxide formation time is of the order of 
milliseconds. These widely different time constants 
present classical methods (such as the popular explicit 
Runge-Kutta method) with the following difficulty: to 
ensure stability of the numerical solution, these methods 
are restricted to using very short time steps that are 
determined by the smallest time constant. However, the 
time for all chemical species to reach near -equilibrium 
values is determined by the largest time constant. As a 


result, classical methods require excessive amounts of 
computer time to solve stiff systems of ordinary differ- 
ential equations (ODE’s). 

Several approaches for the solution of stiff ODE’s 
have been proposed; for details, see the reviews by 
Lomax and Bailey (ref. 3), Seinfeld, et al. (ref. 4), 
Enright and Hull (ref. 5), and Shampine and Gear 
(ref. 6). Of all these techniques the general-purpose codes 
EPISODE and LSODE (refs. 7 to 10) are regarded as the 
best available “packaged” codes for the solution of stiff 
systems of ODE’s. However, although these codes may 
be the best available for solving an arbitrary system of 
ODE’s, it may be possible to construct a superior method 
for solving a particular system of ODE’s governing the 
behavior of a specific problem. In this vein. Young and 
Boris (ref. 11), Pratt (refs. 12 to 15), and Zeleznik and 
McBride (ref. 16) have developed codes for the specific 
purpose of integrating large systems of chemical kinetic 
rate equations. 

The objective of the present investigation is to identify 
the fastest algorithm currently available for the numerical 
integration of combustion kinetic rate equations.! The 
motivation behind this work is the increasing interest in 
(1) modeling the reaction mechanisms describing the 
consumption of fuels and pollutant formation and 
destruction and (2) multidimensional modeling of 
reactive flows, which includes the equations of fluid 
motion. The former results in the need to integrate large 
systems of nonlinear ODE’s (reaction rates). The latter 
results in the need to integrate these rate equations at 
several thousand grid points. To make such calculations 
practicable, it is necessary to have a very fast homo- 
geneous batch chemistry integrator. 

In the present report, currently available techniques are 
examined by application to two test problems drawn 
from combustion kinetics. A detailed comparison of the 
efficiency and accuracy of these techniques is presented, 
and recommendations are made for ways to increase the 
speed of general-purpose codes, as applied to the present 
problem. 


!lti this report, attention is restricted to adiabatic, constant pressure 
(hence, isenthalpic), exothermic chemical reactions. 



Symbols 


Aj, A 


^J’ ^-j 


CPU 

^p,i 

Ei 


Ej, E-j 


EPS 


ERMAX 


fi 


Ho 

HO 

hi 

lERROR 

ITMAX 

ITOL 

J 


preexponential constants in forward and 
reverse rate equations for reaction j (eq. 
(5)), units depend on reaction type 

exponent-on-ten in modified Arrhenius 
preexponential factor, By = logio Aj, 
B_y=logio A_y, arbitrary units 
total CPU time required on IBM 370/3033, s 

constant-pressure specific heat of species i, 
J/kmole K 

estimated local error in independent variable 

/(eq. (17)) 

activation energy in forward and reverse rate 
equations for reaction j (eq. (5)), cal/mole 

for all methods, except EPISODE and 
DASCRU: local relative error tolerance; 
for EPISODE: local relative error tolerance 
for species with initially nonzero mole 
numbers and for the temperature, and local 
absolute error tolerance for species with 
initially zero mole numbers; for DASCRU: 
local relative error tolerance for a variable 
whose magnitude is greater than or equal 
to 10-3, and local absolute error tolerance 
for a variable whose magnitude is less than 
10-3 

relative error tolerance for Newton-Raphson 
iteration for temperature 

rate of formation of species i (eq. (2)), 
kmole i/kg mixture s 

1 -atmosphere molal-specific Gibbs function 
of species i, J/kmole 

initial mass-specific enthalpy of mixture, 
J/kg 

initial step length to be attempted by 
integrator, s 

molal-specific enthalpy of species i, J/kmole 
error control indicator for EPISODE 

maximum number of corrector iterations to 
be attempted by CHEMEQ and CREKID 

error control indicator for LSODE 
Jacobian matrix: for temperature method A 
of size NSxNS: J\\ = dfj/dnky /,A:=1,NS 
(eq. (11)); for temperature method B, of 
size (NS+ l)x(NS-l- 1): J\y^ = dfj/dnk, i,k = 
1,NS (eq. (11)); Ji,t^s + l = W9T, /=1,NS 
(eq. (13)); .^ns + = (dT/dt), 

A: = 1,NS (eq. (14)); /ns + i,ns+i 
idT/dt) (eq. (15)) 


JJ 

^J’ ^-j 
MF 

NFE 

NJE 

NRRC 

NS 

NSTEP 

«/• 

nj 

Arij 

P 

R 

Rj, R -j 
T 

Ep T_j 
AT 


t 

^switch 


Yi 

^abs 

^abs,/ 


^rel 


total number of distinct elementary reactions 
in reaction mechanism 

forward and reverse rate constants for reac- 
tion j (eq. (5)), units depend on reaction type 

integration method to be used by EPISODE 
and LSODE 

temperature exponent in forward and reverse 
rate constants for reaction j (eq. (5)) 

total number of functional (i.e., derivative) 
evaluations 

total number of Jacobian matrix evaluations 

total number of reaction rate constant 
evaluations 

total number of distinct chemical species in 
gas mixture 

total number of steps required to solve 
problem 

mole number of species i, kmole species i/ 
kg mixture 

molecularities of forward and reverse 
reactions j (eq. (12)) 

reciprocal of mixture mean molar mass 
(eq. (10)), kmole/kg 

difference in molecularities of reverse and 
forward reactions j (eq. (6c)) 
pressure, N/m2 

universal gas constant, 8314.3 J/kmole K 
(1.9872 cal/mole K) 

molar forward and reverse rates per unit 
volume for reaction j , kmole/m^ s 

temperature, K 

activation temperatures in forward and 
reverse rate constants for reaction j, 
Tj--Ej/R\ T^j = E-j/R, K 

maximum temperature change allowed before 
reaction rate constants are updated. For 
CREKID, CHEMEQ-B, and Di^SCRU-B 
thermodynamic data are also updated only 
for temperature changes greater than AT, K 
time, s 

time at which error control performed by 
EPISODE is switched from semirelative to 
pure relative, s 

mole fraction of species i 
local absolute error tolerance (eq. (18)) 
local absolute error tolerance for species i (eq. 
(19)) 

local relative error tolerance (eq. (17)) 
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Stoichiometric coefficients for species i in 
forward and reverse reactions j (eq. (1)); 
number of kilomoles i in elementary reaction 
j as a reactant and as a product, respectively 

p mixture mass density, kg/m^ 


Problem Statement 

The (initial value) problem may be stated as follows: 
given (1) a set of initial conditions (/?,• (/= 1,NS), where «, 
is the mole number of species i (kmole species i/kg 
mixture) and NS is the total number of distinct chemical 
species involved in the combustion reaction and the 
temperature T (K)) at time f = 0, (2) the pressure p, and 
(3) the reaction mechanism, find the mixture composition 
and temperature at the end of a prescribed time interval. 

All chemical reactions are assumed to be elementary 
reactions of the type 


/_F \ NS 

=^y7^^xp^— JI 

NS 

= kj'^^(priky^i (4a) 

and 

R^j=A _y7^-^xp ^ j n {pnk)'^\ 

/ — T* A NS 

=A _y7^-^xp Jlj ipnk)'^> 

NS 


NS NS 

D E "ij^/ y=i.JJ (1) 

I = 1 / = 1 

where v{y and vij are the stoichiometric coefficients of 
species i (with chemical formula X,) in reaction j as a 
reactant and as a product, respectively, and JJ is the total 
number of distinct elementary reactions in the given 
reaction mechanism. 

The ordinary differential equations describing 
adiabatic, homogeneous, gas-phase combustion reactions 
are as follows: 


-^=Mrik,T) i,k=UNS \ 

«/(^ = 0) = given ( (2) 

T (/ = 0) = given / 


where 

JJ 

y) = - P - 1 E ("ij - "ijX^y - R -j) (3) 

J=l 

where the molar reaction rates per unit volume Rj and 
R -j are given by 


In equations (2) to (4), p is the mixture density (kg/m3) 
and Ajf A_j, Nj, N_j, Ej, Tj {-Ej/R) E^j, and T-j 
{=E^j/R) are constants in the modified Arrhenius rate 
expressions of the forward Rj and reverse R^j rates of 
reaction j. The forward kj and reverse k^j rate constants 
for reaction j are given by 

kj = A y 7^>exp ^ ^ (5a) 

k^j=A _y T'^-Jexp ^ ^ (5b) 


In this study the reverse rate constants are calculated 
from the forward rate constants and the concentration 
equilibrium constant by using the principle of detailed 
balancing (ref. 17). The resulting equations are 


NS 

E (''ij - 

T-j=Tj+— (6a) 

R 

[ ,NS 

ili J (6b) 
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where hj is the molal-specific enthalpy of species i 
(J/kmole), is the 1 -atmosphere molal-specific Gibbs 
function of species i (J/kmole) and A/iy is given by 

NS 

A«y= (6c) 

;' = 1 

Equality of the temperature exponents in equations 
(5b) and (6b) for k-j together with equation (5a) results 
in the following relation for N^f. 

(6d) 


Some of the techniques used in the present study 
require the evaluation of the Jacobian matrix J 
(</ik = 3/,/ dn^; i,k = 1 ,NS). Differentiation of equation (3) 
with respect to gives the following expression for /jk; 

jj f 

Ak=- iprtk ) ~ ‘ E (''ij “ »'ij) - *'kj-^ -j) + ^ 
y=l 

JJ 

+ (p«w)~^ E (11) 

y = i 

where nj and rtj'are, respectively, the molecularities of the 
forward and reverse reaction j and are given by 


For a constant-pressure, adiabatic combustion reac- 
tion, the following enthalpy conservation equation 
constitutes an algebraic constraint on equations (2) to (4): 

NS 

= = constant (7) 

/ = 1 


NS 

n'j= Y, "ij 
/= 1 


NS 

E "ij 

/= 1 


(12) 


where Hq is the mass-specific enthalpy of the mixture 
(J/kg). 

Time differentiation of equation (7) provides the 
following equation for the time rate of change of 
temperature: 


dt 


NS 

- E 


NS 

E 


( 8 ) 


When temperature is included as an explicit 
independent variable (see the section Evaluation of 
Temperature for more details), the Jacobian elements 
bfi/dT (/ = 1,NS), d/dnu idT/dt) (A:=1,NS), and 
d/bT{dT/dt) are also needed. These are obtained by 
differentiating equation (3) with respect to T and 
equation (8) with respect to and T. These operations 
yield 


dT 





where Cpj is the constant-pressure specific heat of 
species i (J/kmole K). 

The density p of the mixture is given by the equation of 
state for an ideal gas mixture 


P = 


P 

RTtim 


(9) 


-R-,{N.j+ 




NS 

/?1 





NS 

E ^fp,i 


(13) 


(14) 


where p is the absolute pressure (N/m2), R the universal 
gas constant (J/kmole K), and the reciprocal of the mean 
molar mass of the mixture 

NS 

E (10) 

/= 1 


dT\dt) 




NS 

E niCp,i 

l=\ 

(15) 
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Methods Studied 

The methods examined in this study include the 
general-purpose packages EPISODE and LSODE (refs. 7 
to 10), developed for an arbitrary system of ODE’s, and 
the specialized codes CHEMEQ (ref. 11), CREKID (ref. 
15), and GCKP84 (refs. 16 and 18), developed 
specifically to integrate chemical kinetic rate equations. 
In addition, an explicit Runge-Kutta-Merson differential 
equation solver (ref. 19) (IMSL routine DASCRU) is 
used to illustrate the difficulty associated with integrating 
chemical kinetic rate equations by classical methods. 
These methods are summarized in table I and discussed in 
more detail in appendix A. 

The packages EPISODE and LSODE, based on the 
methods of Gear (refs. 20 and 21), consist of a variable- 
step, variable-order implicit Adams method (suitable for 
nonstiff problems) and a variable-step, variable-order 
backward differentiation method (suitable for stiff 
problems). A range of corrector iteration methods— from 
functional iteration to chord method (a variant of 
Newton’s method) with a banded Jacobian generated 
internally— is included in these packages. The user selects 
both the basic method and the corrector iteration method 
by means of a method flag MF. 

In CHEMEQ the ODE’s are separated into two classes, 
stiff and normal, at the beginning of each time step. A 
classical second-order predictor-corrector method is used 
for equations classified as normal. For stiff equations a 
simple asymptotic integration formula is used. 

CREKID is based on the exponentially fitted 
trapezoidal rule developed by Liniger and Willoughby 
(ref. 22) and by Brandon (refs. 23 and 24). The algorithm 
includes special treatment of ill-posed initial conditions 
and automatic selection of functional iteration or 
Newton iteration. 

GCKP84, developed by Bittker and Scullins (ref. 18), is 
a new general chemical kinetics program that supersedes 


their previous GCKP codes (ref. 25). The new code uses 
the integration technique developed by Zeleznik and 
McBride (ref. 16) specifically to integrate chemical 
kinetic rate equations. Details of this integration 
technique are not yet available. 

DASCRU is an explicit fourth-order Runge-Kutta- 
Merson ODE solver. This method requires five derivative 
evaluations per step. The additional derivative evaluation 
provides an estimate of the local truncation error 
(ref. 19). 


Test Problems 

The algorithms summarized in the preceding section 
were applied to two test problems drawn from com- 
bustion kinetics. Both problems include all three regions 
of interest to a combustion researcher: induction, heat 
release, and equilibration. 

Test problem 1, taken from Pratt (ref. 12), describes 
the ignition and subsequent combustion of a mixture of 
33 percent carbon monoxide and 67 percent hydrogen 
with 100 percent theoretical air, at a pressure of 10 
atmospheres and an initial temperature of 1000 K. It 
comprises 12 reactions (table II) that describe the time 
evolution of 11 species. Test problem 2, taken from 
Bittker and Scullins (ref. 18), describes the ignition and 
subsequent combustion of a stoichiometric mixture of 
hydrogen and air, at a pressure of 2 atmospheres and an 
initial temperature of 1500 K. It involves 30 reactions 
(table III) that describe the time evolution of 15 species. 

The initial values for the species mole fractions and 
temperature are given in tables IV and V for test 
problems 1 and 2, respectively. Also given in these tables 
are the equilibrium values for the species mole fractions 
and temperature calculated by using a Gibbs function 
minimization routine (ref. 26). Both test problems were 


TABLE I.— SUMMARY OF METHODS STUDIED 


Code or 
method 

Description 

GCKP84 

Details not yet available 

CREKID 

Variable-step, predictor-corrector method based on an exponentially fitted trapezoidal 
rule; includes filtering of ill-posed initial conditions and automatic selection of 
functional iteration or Newton iteration 

LSODE; 

Variable-step, variable-order backward differentiation method with a generalized 

EPISODE 

Newton iteration^ 

CHEMEQ 

Variable-step, second-order predictor-corrector method with an asymptotic integration 
formula for stiff equations 

DASCRU 

Variable-step, fourth-order, explicit Runge-Kutta-Merson solver 


® Other options are included in these packages. 
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TABLE II.— REATION MECHANISMS AND RATE 
CONSTANTS USED IN TEST PROBLEM 1 


Reaction 

Rate constants^ 

B 

N 

Ey kcal/mole 

C0 + 0H = C02 + H 

11.49 

C 


0.596 

H+02 = 0+0H 

14.34 



16.492 

H2 + 0 = H + 0H 

13.48 



9.339 

H20 + 0 = 0H+0H 

13.92 



18.121 

H+H20 = H2 + 0H 

14.0 

\ 


19.870 

N + 02 = N0 + 0 

9.81 

1.0 

6.250 

N2 + 0 = N + N0 

13.85 

0 


75.506 

N0 + M = N + 0 + M 

20.60 

-1.5 

149.025 

H + H+M-H 2 +M 

18.00 

-1.0 

0 

0 + 0 + M = 02 + M 

18.14 

-1.0 

.34 

H+0H + M = H20 + M 

23.88 

-2.6 

0 

H2 + 02 = 0H + 0H 

13.00 

0 

_ 

43.0 


®Rate constant k= 10^7^cxp(-£//?7). 


TABLE III.— REACTION MECHANISMS AND RATE 
CONSTANTS USED FOR TEST PROBLEM 2 


Reaction 

Rate constants® 

B 

N 

Ey kcal/mole 

H + 02 = 0H + 0 

14.342 

0 

16.790 

0 + H2 = 0H + H 

10.255 

1.0 

8.900 

H2 + 0H = H20 + H 

13.716 

0 


6.500 

0H + 0H = 0 + H20 

12.799 



1.093 

H + O 2 + M = HO 2 + M 

15.176 



-1.000 

0 + 0 + M = 02 + M 

13.756 

f 

-1.788 

H + H + M = H2 + M 

17.919 

-1.0 

0 

H + 0H + M = H20 + M 

21.924 

-2.0 

0 

H2 + H02 = H20 + 0H 

11.857 

0 

18.700 

H202 + M = 0H + 0H + M 

17.068 



45.500 

H2 + 02 = 0H + 0H 

13.000 



43.000 

H + H02 = 0H + 0H 

14.398 



1.900 

0 + H02 = 0H + 02 

13.699 



1.000 

OH+H 02 -H 20+02 

13.699 



1.000 

H02 + H02 = H202 + 02 

12.255 



0 

0H + H202 = H20 + H02 

13.000 



1,800 

0 + H202 = 0H + H02 

13.903 



1,000 

H+H 202 -H 20 +OH 

14.505 



9.000 

H02 + N0 = N02 + 0H 

13.079 



2.380 

0 + N02 = N0 + 02 

13.000 



.596 

N0 + 0 + M = N02 + M 

15.750 



-1.160 

N02 + H=N0 + 0H 

14.462 

\ 


.795 

N + 02 = N0 + 0 

9.806 

1.0 

6.250 

0 + N2 = N0 + N 

14.255 

0 

76.250 

N+OH=NO+H 

13.602 



0 

N20 + M = N2 + 0 + M 

14.152 



51.280 

0-|-N20=^N2 + 02 

13.794 



24.520 

0 + N20 = N0 + N0 

13.491 



21.800 

N+N02 = N0 + N0 

12.556 



0 

0H + N2 = N20 + H 

12.505 



80.280 


®Rate constant /r = 10^7^exp(-£//?7). 


integrated over a time of 1 ms to obtain near 
equilibration of all species. 


Discussion of Results 

In this section, we discuss the computational work 
(which we shall take as a measure of the efficiency of the 
algorithm) and the accuracy of the techniques selected for 
our study. All of the codes were applied to the two test 
problems discussed above. All codes were run on the 
NASA Lewis Research Center’s IBM 370/3033 computer 
using single-precision accuracy, except GCKP84, which 
was in double precision. 

Computational Procedure 

A typical computational run consisted of initializing 
the time (t, set equal to 0), species mole numbers, 
temperature, and the CPU time.2 The integrator was then 
called with values for the necessary input parameters, 
which are discussed in the section Efficiency 
Comparisons, and the elapsed time (1 ms for both 
problems) at which the solution is to be terminated. The 
integrator returns to the main calling program with the 
computed solutions for the mole numbers and, for some 
methods, the temperature. The integrator also returns 
with values for the following parameters, which give a 
measure of the computational work required to solve the 
problem: total number of steps NSTEP, total number of 
functional (i.e., derivative) evaluations NFE, total 
number of Jacobian matrix evaluations (NJE = 0 for 
CHEMEQ and DASCRU), and, for reasons presented 
later, total number of rate constant evaluations NRRC. 
On return from the integrator the computer time CPU 
required to solve the problem was calculated. 

Evaluation of Temperature 

Of the codes tested, only CREKID and GCKP84 were 
written explicitly for the integration of exothermic, non- 
isothermal, combustion rate equations. These therefore 
have built-in procedures for calculating the temperature. 
For the other codes the temperature was computed by 
using one of two different methods, labeled as methods A 
and B.3 


^Before this was done, various preprocessors were called to read in 
thermochemical and reaction rate data and to compute the initial 
mixture properties and the equilibrium composition and temperature. 
This does not affect the work required by the integrator. The storage 
and work requirements of these preprocessors are therefore not 
included in this study. 

^The following convention was adopted in naming these other codes: 
those using temperature method A were given the suffix A (e.g., 
LSODE-A, EPISODE-A, etc.) and those using temperature method B 
were given the suffix B (e.g., CHEMEQ-B, DASCRU-B, etc.). 
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TABLE IV.— COMPOSITIONS AND TEMPERATURES FOR TEST PROBLEM 1 


[Solution generated with LSODE-B and EPS =10 ^.] 


Compo- 

Time, r, s 

nent and 




— 





temper- 

0 

9xl0-« 

10-5 

5x10-5 

10 -“ 

5xl0-“ 

10-3 

00 

ature 









Composition, species mole fraction 

CO 

8.319x10-2 

8.311x10-2 

8.292x10-2 

3.123x10-2 

2.562x10-2 

1.802x10-2 

1.754x10-2 

1.798x10-2 

CO2 

0 

8.359x 10-5 

2.886x10-“ 

6.078x10-2 

6.715x10-2 

7.557x10-2 

7.610x10-2 

7.564x10-2 

H 

0 

1.400x10-2 

4.572x10-2 

5.958x10-2 

3.131x10-2 

1.177x10-2 

1.084x10-2 

1.059x10-2 

H2 

1.664x10-' 

1.641X10-' 

1.588x 10-' 

1.373x10-2 

1.030x10-2 

6.559x10-2 

6.357x10-2 

6.580x10-2 

H2O 

0 

1.519x10-2 

5.221x10-2 

1.608x10-' 

1.682x10-' 

1.762x10-' 

1.767x10-' 

1.768x10-' 

N 

0 

2.419x10-'* 

2.419x10-'* 

3.258x10-2 

4.117x10-7 

3.778x10-2 

3.364x10-2 

2.209x10-2 

NO 

0 

5.501x10-'* 

1.725X10- '2 

6.295x10-5 

1.993x10-4 

9.975x10-“ 

1.666x10-2 

5.347x10-2 

N2 

6.256x 10-' 

6.256x10-' 

6.257X10-' 

6.919x10-' 

6.975x10-' 

7.033x10-' 

7.034x10-' 

7.014x10-' 

0 

0 

1.699x10-“ 

5.678x 10-“ 

3.803x10-2 

2.173x10-2 

8.779X10-“ 

7.986x10-“ 

7.048x10-4 

OH 

0 

4.569x 10-5 

1.508X 10-“ 

1.301x10-2 

1.082x10-2 

7.552x10-2 

7.238x10-2 

6.762x10-2 

O2 

1.248x10-' 

1.239x10-' 

1.217x10-' 

1.878x10-2 

1.484x10-2 

9.752x10-2 

9.171x10-2 

7.826x10-2 

Temperature, T, K 

Ta 

1000 

1001 

1006 

2407 

2512 

2623 

2629 

2619 

Tb 

— 

1001 

1006 

2407 

2512 

2623 

2628 

— 


TABLE V.— COMPOSITIONS AND TEMPERATURES FOR TEST PROBLEM 2 


[Solution generated for LSODE-B and EPS =10 ^.] 


Compo- 

Time, /, s 

nent and 










temper- 

0 

3x10-* 

5x10-6 

10-5 

0 

X 

10'“ 

0 

X 

10-3 

00 

ature 










1 Composition, species mole fraction 

Ar 

6.571 X 10-3 

6.572x 10-3 

6.605 X 10-3 

6.744 X 10-3 

7.080x10-3 

7.204x10-3 

7.323x10-8 

7.324x10-3 

7.326X 10-3 

CO2 

2.110x10-4 

2.111x10-4 

2.121x10-4 

2.166x10-'* 

2.274x10 -4 

2.314x10-4 

2.352x10-4 

2.352x10-4 

2.353x10-“ 

H 

0 

3.006X 10-3 

1.073x 10-1 

9.848x10 2 

4.542x10-2 


1.870x 10-2 


1.853x10-2 

HO2 

0 

3.261 X 10-5 

1.183x10-5 

5.878x 10-6 

5.121x10-6 

8.674x10-6 

1.447x10-5 

1.398x10-5 

1.391x10-5 

H2 

2.951 X 10-1 

2.896 X 10-1 

9.272x 10-2 

6.201x10-2 

6.063x10-2 

5.544x 10-2 

4.924x10-2 

4.966x10-2 

4.974x10-2 

H2O 

0 

3.807 X 10-2 

1.429x10-1 

1.793x10-' 

2.171x10' 

2.353x10-' 

2.549x10' 

2.549x10-' 

2.549x10-' 

H2O2 

0 

3.442x10-8 

3.790x 10-* 

5.603x10-* 

2.629x10'* 


6.377x10-* 

6.185x10-* 

6.160x10-* 

N 

0 

1.718x10-12 

2.039 X 10-10 

4.927X10-* 

3.084x10'* 

9.124x10-^ 


4.067x10-6 

4.144x10-* 

NO 

0 

4.894X 10-1? 

3.769X 10-9 

2.925x 10-'^ 

1.428x10'“ 

8.095x 10-4 

7.615x10-3 

9.752x 10-3 

1.008x10-2 

NO2 

0 

1.639X 10-15 

2.112X 10-1^ 

3.520x 10-11 

1.928x10-8 

1.374x10-* 

1.568x 10-6 

1.973x 10-6 

2.035x10'* 

N2 

5.501x10-1 

5.502 X 10-1 

5.530x 10-1 

5.646x10' 

5.926x10-' 


6.094x10-' 

6.084x10-' 

6.082x10-' 

N2O 

0 

1.090x10-* 

6.095x 10-8 

1.149x10-* 

1.704x10-'^ 

3.010x10-* 

6.382x10-* 


6.345x10-* 

0 

0 

6.325x 10-4 

2.533x 10-2 

2.841x10-2 

1.690x10-2 

1. 204 X 10-2 

7.113x10-3 

6.866x10-3 

6.833x 10-3 

OH 

0 

3.358x 10-4 

1.452x 10-2 

2.440x10-* 

3.495x 10-2 

3.485 X 10-2 

3.068x 10-2 

3.011x10-2 

3.003x10-* 

O2 

1.480x10-1 

1.456x 10-1 

5.738x 10-2 

3.583x10-* 

2.495x 10-2 

2.078x10-* 

1.480x10-2 

1.417x10-* 

1.408x10-* 

Temperature, T, K 


1500 



WtBM 
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In method A the temperature was calculated from the 
mole numbers and the initial mixture enthalpy by using 
the algebraic enthalpy equation (eq. (7)). A Newton- 
Raphson iteration technique (with a user-specified 
relative error tolerance ERMAX) was used to solve 
equation (7) for the temperature. In this method the 
temperature does not enter into the problem as an explicit 
independent variable, so the number of independent 
ode’s is equal to the number of species NS and the 
Jacobian matrix is of size NS by NS. The integrator 
therefore tracks only the solution for the mole numbers. 
As shown later, this did not introduce significant errors 
because the species concentrations changed faster than 
the temperature (figs. 1 and 2) and the step length was 
determined by the rate at which the variables changed. 
Since the integrator did not compute a solution for the 
temperature, it was evaluated from the solution for the 
mole numbers returned by the integrator. In addition, the 
temperature was computed whenever the species time 



10-6 10-5 10-4 10-3 


Time, t, s 

Figure I.— Variation with time of species mole fractions and 
temperature for test problem 1, Solution generated with LSODE-B 
and EPS = 10-5. 


derivatives (dn,/dt) and the Jacobian matrix 
(/jl( = 3/)/3njt,‘ /,A: = 1,NS) were evaluated. 

In method B the temperature was treated as an 
additional independent variable and evaluated by 
integrating its time derivative (eq. (8)). This increased the 
number of independent ODE’s to NS-)-l, and the 
computation of the Jacobian matrix (of size (NS-t- 1) by 
(NS + 1)) involved the calculation of 2NS-I-1 additional 
terms. In this method, the integrator tracks the solutions 
for both the temperature and the species mole numbers. 

Accuracy Estimate 

To compare the accuracy of the algorithms, standard 
solutions must be established since exact solutions for the 
test cases are not known. The solutions used as standards 
were those generated by using LSODE-B (which is 
LSODE using method B for computing the temperature) 
because LSODE and its predecessors have been 



Time, t, s 


Figure 2.— Variation with time of species mole fractions and 
temperature for test problem 2. Solution generated with LSODE-B 
and EPS= 10-5. 
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extensively tested on a wide variety of problems (refs. 7, 
8, 27 to 29). A low relative error tolerance (EPS = 10 "5) 
was used, and the method flag MF was set equal to 21 
(stiff method, user-supplied analytic evaluation of the 
complete Jacobian matrix). As discussed in appendix B, 
LSODE requires the specification of values for both the 
local relative and absolute error tolerances. In the test 
problems examined in this study the species mole 
numbers and temperature differed widely (tables IV and 
V), so relative error control was appropriate. Pure 
relative (or absolute) error control can be realized by 
using a value of zero for the absolute (or relative) error 
tolerance. However, since some of the mole numbers had 
zero initial values, a truly relative error criterion could 
not be used. To make the local error control mostly 
relative, low values for the absolute tolerance for species 
mole numbers and zero for temperature were used. 
Values for the absolute error tolerances for the species 
mole numbers were obtained by progressively decreasing 
them until the temperature-time trace showed essentially 
no change with a further decrease. The values used for 
the absolute error tolerances for the species mole 
numbers were 10- 14 and 10- H for test problems 1 and 2, 
respectively. 

To determine if the use of solutions generated with 
LSODE as standards biases the test in favor of LSODE, 
we explore the use of other standards in the section 
Efficiency Comparisons and show that these do not give 
different results. 

The standard solutions (for the mole fractions and 
temperature) are shown in figures 1 and 2 for test 
problems 1 and 2, respectively. The values for the mole 
fractions and temperatures at various times are shown in 
tables IV and V for problems 1 and 2, respectively. The 
third column in these tables presents immediate 
postinduction (i.e., end of the essentially isothermal 
reaction period) values, and the remaining columns show 
development of postinduction profiles. The last column 
shows the equilibrium solution obtained by using a Gibbs 
function minimization routine (ref. 26). 

Although these standard solutions were generated by 
using temperature method B, method A was also used 
during these computations to evaluate the temperature, 
using only the solution for the mole numbers returned by 
the integrator. This was done to see if any consistent 
differences could be observed between values generated 
by using method A (T^ in tables IV and V) and those 
generated by using method B (T^in tables IV and V). The 
consistently excellent agreement between and Tq 
indicates that either method A or B would suffice for 
computing the temperature. 

Storage Requirement and Startup Time 

Multidimensional models for reactive flows require the 
integration of large sets of reaction rate equations at 


several thousand grid points for relatively short times 
between fluid mechanic time steps (ref. 1 1). These models 
generally also have large storage requirements. Hence 
reaction rate integration techniques with both a low 
storage requirement and a low initialization (startup) 
time are needed. 

The storage requirements (for 20 species and 36 
reactions) for the single-precision versions of CREKID, 
LSODE, EPISODE, CHEMEQ, and DASCRU are given 
in table VI. The storage requirements for LSODE, 
EPISODE, and DASCRU are for versions that included 
both methods A and B for calculating the temperature. 
The temperature method to be used was specified via a 
temperature-method flag. Also, for LSODE the given 
figure does not include the storage required by 
subroutines (included in this package) that are not 
essential for its execution with method flag MF = 21. The 
storage shown in table VI for GCKP84 is that required by 
the double-precision version. We note that the special- 
purpose codes CHEMEQ and CREKID required much 
less storage than the general-purpose codes EPISODE 
and LSODE. GCKP84 required more storage than the 
other codes because of the precision used and also 
because of the various options built into it. Although 
GCKP84 is a special-purpose code in that it has been 
developed for chemical kinetics problems, it can be used 
to solve a wider variety of problems than CHEMEQ and 
CREKID (ref. 18). 

The CPU time required by each code to successfully 
complete the first step was taken as a measure of the 
startup time. Each code was run with a low value for the 
elapsed time to ensure that only one step was taken to 
complete the problem. Note, however, that because of 
the automatic filtering of initial conditions (appendix A), 
CREKID took two steps to return with a solution. The 
CPU times (in milliseconds) required for the first step 
(and the first two steps for CREKID) are given in table 
VII. The general-purpose codes required longer 
initialization times than the special-purpose codes. 
GCKP84 required a much longer startup time than the 
other codes for the same reasons that its storage 
requirement is greater and also because of the extra work 


TABLE VI.— STORAGE 
REQUIREMENTS 


Method 

Storage size, 
bytes 

GKCP84® 

205 948 

CREKID 

26 304 

LSODE 

47 152 

EPISODE 

35 116 

CHEMEQ-A 

10 832 

CHEMEQ-B 

10 552 

DASCRU 

10 928 


^Double-precision version. 
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TABLE VII.— CENTRAL PROCESS- 
ING UNIT TIMES REQUIRED TO 
SUCCESSFULLY COMPLETE 
FIRST STEP 


Method 

Problem 1 

Problem 2 

CPU time, CPU ms, 

GCKP84 

32 

44 

CREKID 

a9 

HI 

LSODE-A 

20 

24 

LSODE-B 

19 

24 

EPISODE-A 

14 

18 

EPISODE-B 

13 

18 

CHEMEQ-A 

6 

10 

CHEMEQ-B 

6 

8 

DASCRU-A 

7 

11 

DASCRU-B 

6 

10 


®Time required for first two steps. 


involved in handling species whose mole numbers have 
zero initial values (see ref. 18 for details). 

Efficiency Comparisons 

The procedure described in the section Computational 
Procedure was used to study the computational work 
required by the different techniques examined in this 
study. For LSODE, EPISODE, CHEMEQ, and 
DASCRU both temperature methods A and B were 
attempted. 

The codes examined in the present study require the 
specification of several parameters in addition to the 
local error tolerance EPS required of the solution. Values 
for these parameters that minimized the computational 
work required by each code were obtained by using a 
trial-and-error process. Following are the user -supplied 
parameters (excluding the local error tolerance EPS and 
the elapsed time at which the integration is to be 
terminated) required by each code (see appendix B for a 
detailed discussion) — for EPISODE; the method flag 
MF, the error control to be performed lERROR, and the 
initial step length HO to be used; for LSODE: the method 
flag MF, the error control to be performed ITOL, and 
the values for the absolute error tolerances for the 
independent variables; for CHEMEQ: the maximum 
number of corrector iterations ITMAX allowed before 
nonconvergence is declared; for CREKID: the maximum 
number of corrector iterations ITMAX allowed before 
nonconvergence is declared and the maximum 
temperature change AT allowed before thermodynamic 
properties and reaction rate constants are updated; for 
GCKP84: since details for this technique are not yet 
available, default values for all parameters to be used; for 
DASCRU: a guess for the initial step length HO to be 
used. 

The following selection procedure was adopted in 
using these techniques: each code was run with a value of 


10-2 for the error tolerance EPS, and a solution (with 
output at various times) was generated. The computed 
temperature-time profile was then compared with the 
standard solution and was accepted if within 50 kelvins. 
Otherwise, it was rejected and a lower value for EPS was 
tried. This ensured that computational work was 
compared with codes of comparable accuracy. All of the 
results presented herein, except EPISODE and GCKP84 
for test problem 1, satisfied this criterion. In other words, 
for these codes no EPS > 10-6 resulted in acceptable 
agreement. With DASCRU the temperature did not 
satisfy this criterion during early heat release for test 
problem 1. 

For test problem 2, some runs with DASCRU and 
LSODE predicted zero mole numbers for the nitrogen 
dioxide at times when the standard solutions had risen to 
measurable levels (of the order of 0.1 ppm). All codes 
were therefore required to satisfy the following 
additional criterion: a run was accepted only if it 
predicted nonzero mole fractions for all species whose 
standard solution values were greater than or equal to 
10-2. All results presented herein satisfied this criterion. 

Each code was run with the maximum (and lower) 
values of EPS that resulted in acceptable agreement with 
the standard solution. The computational work was 
obtained by following the procedure outlined earlier. 
Figures 3 and 4 present the computational work 
(expressed as the CPU time in seconds required on the 
NASA Lewis Research Center’s IBM 370/3033 com- 
puter) plotted against the local error tolerance EPS for 
test problems 1 and 2, respectively. For all codes except 
EPISODE and DASCRU, EPS is the local relative error 
tolerance required of the numerical solution. For 
EPISODE, EPS is a mixed relative and absolute error 
tolerance— relative for species with initially nonzero mole 
numbers and for temperature (method B) and absolute 
for species with initially zero mole numbers'*. For 
DASCRU, EPS is absolute for variables whose 
magnitude is less than 10-3 and relative otherwise. For 
this study, because the maximum permissible local 
relative error in the temperature calculated by using 
method B was EPS, ERMAX (the relative error allowed 
in the Newton-Raphson iteration procedure used in 


^Because some mole numbers had zero initial values, pure relative 
error control (option IERROR = 2) could not be used with EPISODE. 
The option IERROR = 3 was used, instead. This is a semirelative error 
control where, for a variable that is initially nonzero, the error control is 
relative. For a variable that is initially zero, the error control is absolute 
until the variable reaches 1 in magnitude, when the control becomes 
relative. Since none of the mole numbers reaches a value of unity, the 
error control is always absolute for species with initially zero mole 
numbers. In CHEMEQ and CREKID, mole numbers less than 10"^® 
were set equal to 10“^^. In addition, in CREKID the convergence test 
was applied only to species whose mole numbers are greater than 10 “ 20 
(appendix A). 
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LSODE-B 
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EPISODE-B 
CHEMEQ-A 
CHEMEQ-B 
DASCRU-A 
DASCRU-B 




10 '^ 10"5 10'3 10 “^ 

Error tolerance. EPS 

Figure 3.— Variation of the CPU time with error tolerance for test problem 1, All runs on IBM 370/3033, 


— V— GCKP84 
— O— CREKID 
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■ LSODE-B 
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■■ EPISODE-B 
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dascru-a 

— ^1— DASCRU-B 



□ 


Error tolerance, EPS 

Figure 4. — Variation of the CPU time with error tolerance for test problem 2. All runs on IBM 370/3033. 



method A to solve the algebraic energy equation) was set 
equal to EPS to make the two methods comparable. For 
the same reason, with LSODE-B the absolute error 
tolerance for the temperature was set equal to zero. 

For test problem 1, very small values for the error 
tolerance had to be used for GCKP84, EPISODE, and 
DASCRU (fig. 3). Even with EPS = 10-6, EPISODE-A, 
EPISODE-B, and GCKP84 did not accurately track the 
temperature during ignition and heat release, when the 
solution changed rapidly (fig. 1). For EPS > 5x 10-6, 
EPISODE predicted physically meaningless results — little 
or no change from initial values after an elapsed time of 
1 ms. However, with GCKP84 and DASCRU, although 
higher values of EPS resulted in poor solutions during 
ignition and heat release, the predicted solutions during 
equilibration were satisfactory. For example, with 
GCKP84 and DASCRU, values of EPS as high as 
5 X 10-3 and 10-3, respectively, were adequate to track 
the temperature in this regime with an error of less than 
1 percent. The run with GCKP84 and EPS = 10-2 
exhibited serious instability and so was terminated. 

Similar remarks apply to test problem 2 where, again, 
small values of EPS had to be used with EPISODE-A, 
EPISODE-B, and GCKP84 to track correctly the 
temperature during ignition and heat release. During 
equilibration, higher values for EPS could be used 
without incurring error penalties. With EPISODE and 
GCKP84, values of EPS equal to 10-4 and 10-2, 
respectively, were adequate to follow the temperature 
within 1 percent in this regime, although these runs did 
not satisfy the accuracy criterion during ignition and heat 
release. EPISODE-A and EPISODE-B predicted little or 
no change in the composition and temperature after an 
elapsed time of 1 ms for EPS greater than or equal to 
5x10-4 and 5x10-3, respectively. Although the runs 
with EPISODE-B and EPS of 5x 10-4 and 10-3 were 
successfully completed, the solutions were considerably 
inaccurate. Also, these runs were less efficient than the 
run with EPS = 5x10-5. For a more detailed study on 
the variation of the CPU time with EPS, see references 30 
and 31. 

The selection procedure used above required that alt 
integrators track the solutions generated with LSODE. 
This accuracy criterion is subject to the following two 
objections: (1) it biases the test in favor of LSODE at the 
expense of the other integrators, and (2) there is no 
guarantee that LSODE is more aceurate than the other 
codes examined in this study. The other integrators were 
therefore being required to track not the true solutions 
(which are not known) but solutions that could be less 
accurate than what they themselves generated. In other 
words, the use of solutions generated with another 
integrator might have resulted in the acceptance of some 
runs that were rejected above. This possibility was 
explored by further testing. 


For codes5 that failed to satisfy the comparison with 
LSODE, new standards were established. For every one 
of these codes the solutions used as the new standards 
were those generated by the code itself using a low value 
for EPS. All runs that had previously been rejected were 
then compared with these new standards, which were all 
found to be essentially the same as the standard solutions 
generated with LSODE. More specifically, the following 
comparisons were made: for test problem 1 the runs with 
GCKP84 and EPS > 5 X 10-6 were compared with the 
solution generated with GCKP84 and EPS = 10-8, and 
the runs with CHEMEQ-B and EPS > 5x10-3 were 
compared with the solution generated with CHEMEQ-B 
and EPS =10-5. For test problem 2 the runs with 
GCKP84 and EPS s 5 x 10-5 were compared with the 
solution generated with GCKP84 and EPS = 10-8, and 
the runs with EPISODE-A and EPISODE-B with 
EPS ^ 5 X 10“5 were compared with solutions generated 
with EPISODE-A and EPISODE-B, respectively, and 
EPS = 10-6. In making these new comparisons the 
selection criterion remained the same: namely, accept the 
run if the predicted temperature-time profile is within 50 
kelvins of the new standard. These additional 
comparisons did not result in the acceptance of any run 
that had previously been rejected. Furthermore, for test 
problem 1 the run with GCKP84 and EPS = 10-6 failed 
the 50 kelvin requirement when compared with GCKP84 
and EPS = 10-8. These additional comparisons not only 
support the use of solutions generated with LSODE as 
standards, but also imply that LSODE, CREK ID, and 
CHEMEQ-A are more accurate than the other codes. 

Figures 3 and 4 illustrate the difficulty associated with 
using a classical method (in this case the explicit Runge- 
Kutta method used in DASCRU) to integrate combustion 
kinetic rate equations. The CPU times required for the 
two test problems were of the order of 1 and 15 min, 
respectively. Using this technique would make 
multidimensional modeling of practical combustion 
devices prohibitively expensive. 

For test problem 1 the difference in computational 
work required by methods A and B was small (fig. 3), 
with method B being more efficient. For test problem 2 
(fig. 4) the difference was small for large values of EPS. 
But for small values of EPS the difference was more 
marked, with method A being significantly superior to 
method B. The temperature-time profile was steeper for 
test problem 1 (figs. 1 and 2), indicating a stronger 


^These inclu<ie GCKP84 and CHEMEQ-B for test problem 1 and 
GCKP84, EPISODE-A, and EPISODE-B for test problem 2. 
Although some runs with DASCRU- A and DASCRU-B also failed to 
satisfy the comparisons with LSODE, the new tests were not applied to 
DASCRU. The main purpose in using DASCRU was to illustrate the 
difficulty associated with using classical methods to integrate 
combustion kinetic rate equations. 
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coupling between the species and the temperature. This 
may explain why the inclusion of the temperature as an 
additional independent variable worked well for test 
problem 1. But for test problem 2 the additional work in 
computing the temperature rate and the temperature- 
dependent terms in the Jacobian matrix did not lead to 
increased efficiency. 

LSODE and CREKID were superior to the other codes 
examined in this study (figs. 3 and 4). Comparing 
EPISODE and GCKP84 with CREKID and LSODE was, 
however, difficult for the reasons discussed below. The 
error control performed by EPISODE is different from 
that performed by the other codes. As shown later, for 
problems of the type examined in this study, the error 
control used in EPISODE is inferior to that used in 
LSODE. To attain comparable accuracy at early times, 
EPISODE must use much lower values of EPS than 
LSODE, a result seen earlier. Nevertheless, when interest 
was restricted to computational efficiency, EPISODE 
was an attractive alternative to LSODE and CREKID, 
especially for test problem 2 (see ref. 30 for details). In 
using EPISODE, however, a word of caution is in order. 
The computational work can depend strongly on the 
value for the initial step length HO chosen by the user. A 
poor guess for HO can make EPISODE prohibitively 
expensive to use, as shown for test problem 2 in table 
VIII. Note an order -of-magnitude increase in the CPU 
time for a change in HO from 10-7 to 10-8. Although 
not shown here, the solution was found to be adversely 
affected by a poor choice for HO. Also, some values for 
HO resulted in problems with solution instability. 

Comparing GCKP84 with CREKID and LSODE was 
difficult because of the much lower values of EPS used 
by GCKP84. These low values were necessary to satisfy 
the two independent accuracy criteria used above. For 
test problem 1, GCKP84 was significantly slower than 
EPISODE for the same value of EPS. For test problem 2 
its speed was comparable to those of LSODE-B and 
EPISODE-B for the same values of EPS. For larger 


TABLE VIII.— EFFECT OF INITIAL STEP LENGTH ON 
WORK REQUIRED BY EPISODE-A (EPS = 10-5) poR 
TEST PROBLEM 2 


Initial 

step 

length, 

HO, 

s 

Total 
number 
of steps 
required, 
NSTEP 

Total 
number 
of functional 
evaluations, 
NFE 

Total 

number of 
Jacobian 
matrix 
evaluations, 
NJE 

Total 

CPU 

time 

required, 

CPU, 

s 

10-5 

129 

237 

33 

0.79 

10-« 

129 

231 

31 

.78 

10-7 

126 

225 

36 

.79 

10-* 

1168 

2355 

353 

7.91 

10-’ 

1170 

2394 

362 

8.04 

10-10 

133 

231 

32 

.77 


values of EPS, however, GCKP84 was significantly 
slower than LSODE, EPISODE, and CREKID (ref. 30). 
It should be emphasized that GCKP84 is a general- 
purpose chemical kinetics code designed to solve a variety 
of chemical kinetics problems. Consequently the 
overhead associated with functional and Jacobian 
evaluations is higher for GCKP84. In spite of the extra 
work required per step GCKP84 has been shown to be an 
efficient code for performing a wide variety of chemical 
kinetics calculations (ref. 18). 

Computational Tactics 

A simple way of increasing the efficiency of the 
algorithms as applied to the present problem was 
explored. Experience with CHEMEQ-A showed that 
computing the rate constants kj [=AjT^j exp (-7}*/ 7)] 
and k^j [=A_jT^-J exp {-T^j/T)] every time the 
species time derivatives are evaluated is quite inefficient. 
Note that evaluating the rate constants necessitates 
computing both the exponential terms in the expressions 
for kj and k-j (eqs. 5(a) and 6(b)) and the fifth-order 
polynomial expression used for the Gibbs functions (ref. 
26). To reduce the computational work associated with 
evaluating the rate constants, these were updated only 
when the temperature change exceeded a value AT. 
Shown in table IX are the CPU times for various values 
of AT using CHEMEQ-A on test problem 1. Note the 
substantial decrease in computational work — a value of 
AT=0.1 kelvin resulted in about a 40-percent decrease in 
both the number of steps and the number of function 
evaluations. Furthermore the number of rate constant 
evaluations NRRC decreased from 27 736 to 2578 for 
AT =0.1 kelvin. These reductions in computational work 
resulted in about a 70-percent decrease in the CPU time 
required to solve the problem. 

TABLE IX.— EFFECT OF MAXIMUM TEMPERATURE 
CHANGE ALLOWED BEFORE REACTION RATE 
CONSTANT UPDATE ON WORK REQUIRED BY 
CHEMEQ-A (EPS = 10-3) FOR 
TEST PROBLEM I 


Maximum 

temperature 

change, 

AT, 

kelvin 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number of 
functional 
evaluations, 
NFE 

Total 

number of 
reaction rate 
constant 
evaluations, 
NRRC 

Total 
CPU time 
required, 
CPU, 
s 

(a) 

13 272 

27 736 

27 736 

28.4 

0 

13 272 

27 736 

18 624 

24.8 

.05 

8 061 

17 131 

3 693 

10.4 

.1 

8 030 

17 103 

2 578 

10.0 

.5 

13 582 

28 394 

953 

14.4 

1 

13 974 

29 300 

627 

14.7 

(b) 

8 027 

17 107 

2 827 

9.3 


®Not used. 

^Calculated by using eq. (16). 
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In selecting a value for AT, care must be taken to avoid 
poor approximations in the resulting reaction rates, 
which leads to excessive computational work. Table X 
presents an example of such a case. A value of AT" as 
small as 1 kelvin resulted in an order-of-magnitude 
increase in the CPU time required to solve the problem. 
The selection of an optimum value for AT — defined as 
that value which results in minimum computational 
work — is therefore a trial-and-error process, with it being 
a function of the problem and the error tolerance 
specified. 

To avoid repeated runs of the program in search of the 
optimum value for AT, an approximation for it was 
developed. The results are given below and details are 
presented in appendix C. By requiring that the maximum 
relative error in the resultant reaction rates not exceed the 
required relative error tolerance EPS, the following 
expression for AT was derived: 


EPST 


AT= 


max 

j 


+AT' 
T J’ 



( 16 ) 


where T is the current value of the temperature, the bars 
I I denote absolute value, and the maximum is taken over 
all forward and reverse reactions. 

Every time the reaction rate constants were evaluated, 
which occurred only when the temperature change (since 
the last update of the rate constants) exceeded AT, a new 
value of r could be calculated from equation (16). Thus 
equation (16) provides a simple expression for the 
automatic evaluation (through the history of the 
problem) of AT. The computational work resulting from 


the use of equation (16) is given in tables IX and X. 
Although using equation (16) did not necessarily result in 
the fastest algorithm, it did result in shorter CPU times. 
For CHEMEQ-A this decrease was significant; but for 
reasons discussed below, it was small for LSODE-A. 

Using equation (16) resulted in significant savings for 
CHEMEQ and DASCRU for all values of EPS used in 
this study. But with EPISODE and LSODE the savings 
were found to be small, especially when low values of 
EPS were used. The decreases in CPU time were more 
significant for CHEMEQ and DASCRU because these 
two codes require many more reaction rate evaluations 
(NRRC in tables XI and XII) than EPISODE and 
LSODE. For EPISODE and LSODE the computational 
work required for evaluating AT by equation (16) is 
therefore a greater fraction of the work saved by not 
updating the rate constants than for CHEMEQ and 
DASCRU. In addition, when EPS was small, the 
resulting AT has such small values that the decrease in 
NRRC (and hence in computational work) is not 
sufficient to offset the work required to update AT. This 
update occurs more often as EPS is decreased. 

In incorporating equation (16) into codes that use 
temperature method B, the following attempt at further 
reducing the computational work was made. In addition 
to updating the rate constants only for temperature 
changes greater than AT, the thermodynamic properties 
hi and ,• (eq. (8)) were updated only for temperature 
changes greater than AT. This reduced the work 
associated with computing the fifth-order polynomial 
approximations used in evaluating hj and Cpj (ref. 26). 
For CHEMEQ-B and DASCRU-B this resulted in 
increased efficiency for all values of EPS used in this 
study. But for LSODE-B and EPISODE-B no consistent 
differences could be observed. Hence with LSODE-B 
and EPISODE-B the thermodynamic properties were 


TABLE X.— EFFECT OF MAXIMUM TEMPERATURE 
CHANGE ALLOWED BEFORE REACTION RATE 
CONSTANT UPDATE ON WORK REQUIRED 
BY LSODE-A (EPS= lO-"*) FOR 
TEST PROBLEM I 


Maximum 

temperature 

change, 

Ar, 

kelvin 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number of 
functional 
evaluations, 
NFE 

Total 

number of 
Jacobian 
matrix 
evaluations, 
NJE 

Total number 
of reaction 
rate constant 
evaluations, 
NRRC 

Total CPU 
time 

required, 

CPU, 

s 

(a) 

206 

293 

33 

326 

0.63 

0 

206 

293 

33 

263 

,58 

.05 

227 

317 

32 

203 

,58 

,1 

220 

296 

34 

171 

.54 

.5 

438 

642 

109 

209 

1.12 

1 

1224 

2163 

550 

317 

3.88 

(b) 

213 

289 

31 

224 

, *57 


®Not used. 

^’Calculated by using eq. (16). 
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TABLE XL— WORK REQUIRED FOR TEST PROBLEM 1 


[Given CPU times represent minimum time required by each code to solve test problem 1.) 


Method 

Error 

tolerance, 

EPS 

Relative 

error 

tolerance, 

ERMAX 

Initial 

step 

length, 

HO, 

s 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number of 
functional 
evaluations, 
NFE 

Total number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total number 
of reaction 
rate constant 
evaluations, 
NRRC 

Total CPU 
time 

required, 

CPU, 

s 

GCKP84 

10-« 

(a) 

bio-6 

309 

684 

54 

684 

3.13 

CREKID 

10-2 

(a) 

(a) 

84 

280 

32 

91 

.23 

LSODE-A 

10-2 

10-2 

(a) 

114 

165 

27 

105 

.32 

LSODE-B 

5x10-3 

CO 

(a) 

101 

156 

27 

113 

.32 

EPISODE-A 

10-6 

10-2 

10-10 

244 

463 

41 

435 

.70 

EPISODE-B 

10-6 

(a) 

10-’ 

248 

460 

36 

447 

.68 

CHEMEQ-A 

10-2 

10-2 

(a) 

6 741 

13 963 

0 

735 

6.58 

CHEMEQ-B 

10-3 

(a) 

(a) 

11 884 

24 943 

1 

2 638 

9.82 

DASCRU-A 

10-5 

10-2 

10-2 

12 997 

70 075 

1 

31 395 

43.2 

DASCRU-B 

10-5 

(a) 

10-2 

12 822 

69 185 

t 

31 827 

35.5 


®Not needed, 

^Default value. 

Absolute error tolerance for temperature. 


TABLE XIL— WORK REQUIRED FOR TEST PROBLEM 2 


[Given CPU times represent minimum time required by each code to solve test problem 2.] 


Method 

Error 

tolerance, 

EPS 

Relative 

error 

tolerance, 

ERMAX 

Initial 

step 

length, 

HO, 

s 

Total 
number 
of steps 
required, 
NSTEP 

Total 
number of 
functional 
evaluations, 
NFE 

Total number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total number 
of reaction 
rate constant 
evaluations, 
NRRC 

Total CPU 
time 

required, 

CPU, 

s 

GCKP84 

10-5 

(a) 


137 

313 

34 

313 

3.06 

CREKID 

10-3 

(a) 

— 

140 

439 

138 

214 

1.04 

LSODE-A 

10-2 

10-2 

KH 

87 

137 

26 

87 

.47 

LSODE-B 

10-2 

CO 

(a) 

72 

126 

22 

75 

.45 

EPISODE-A 

10-5 

lO-"* 

io-« 

127 

227 

31 

229 

.71 

EPISODE-B 

10-5 

(a) 

io-» 

145 

303 

34 

273 

.86 

CHEMEQ-A 

10-2 

10-2 

(a) 

8 167 

17 033 

0 

1 171 

13.7 

CHEMEQ-B 

10-2 

(a) 

(a) 

8 745 

18 213 


974 

12.0 

DASCRU-A 

10-5 

10-2 

io-» 

94 596 

594 640 


11 666 

400 

DASCRU-B 

10-5 

(a) 

io-» 

94 627 

593 940 


11 057 

310 


^Not needed, 

^Default value. 

Absolute error tolerance for temperature. 


updated every time the derivatives and the Jacobian 
matrix were evaluated. But with CHEMEQ-B and 
DASCRU-B these properties were evaluated only when 
the rate constants were updated. 6 We note here that 
CREKID allows for a user-specified Ar and both the rate 
constants and the thermodynamic properties are updated 
only for temperature changes greater than Ar. 

The changes discussed above were incorporated into 
CHEMEQ, DASCRU, EPISODE, and LSODE. The 
computational work required by each code for test 


6with method A, since a Newton-Raphson iteration technique was 
used to compute the temperature, the thermodynamic properties were 
updated every time the derivatives and the Jacobian matrix were 
evaluated. 


problems 1 and 2 is presented in tables XI and XII, 
respectively. For temperature method A, ERMAX was 
allowed to be different from EPS (cf. figs. 3 and 4) and 
the value resulting in the least computational work was 
used. For LSODE-B the value for the absolute error 
tolerance for the temperature was chosen similarly. 
Tables XI and XII present the values for the user- 
specified parameters that result in the least 
computational work. The CPU times given in tables XI 
and XII represent the minimum time required by each 
code to solve test problems 1 and 2, respectively. 

Comparing tables XI and XII with figures 3 and 4 
shows the savings realized by using equation (16). For 
DASCRU (classical Runge-Kutta method) and 
CHEMEQ the decreases in CPU time were significant. 
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ranging from about 30 to about 65 percent. For LSODE 
and EPISODE the decreases were more modest, ranging 
from negligibly small to about 20 percent. 

Tables XI and XII show that CREKID was the fastest 
code for test problem 1, and LSODE for test problem 2. 
Note that although CREKID is only marginally faster 
than LSODE for problem 1 , it is significantly slower for 
problem 2. We therefore conclude that LSODE is the 
fastest code currently available for integrating chemical 
kinetic rate equations. LSODE is faster than CREKID 
because of (1) better step length control, especially at 
early times (see the section Step Length Comparisons for 
details) and (2) less frequent update of the Jacobian. For 
test problem 2, CREKID requires 138 Jacobian 
evaluations, whereas LSODE-B requires only 22 
Jacobian evaluations (table XII). The main difficulty in 
using LSODE is the trial-and-error procedure necessary 
to obtain optimum values for the absolute error 
tolerances for a given value of EPS. 


Accuracy Comparisons 

The standard solutions established in the section 
Accuracy Estimate and given in tables IV and V were 
used to compare the accuracy of the algorithms examined 
in this study. Tables XIII and XIV give the differences 
from the standard solutions, in percent, for the solutions 
obtained with the optimized codes. The values chosen for 
the user-specified parameters are given in tables XI and 
XII for test problems 1 and 2, respectively. As expected, 
all temperature predictions were in close agreement with 
the standard solutions, except for the immediate 
postinduction solutions with GCKP84 and EPISODE. 
Unlike the other codes CHEMEQ tracked the 
temperature more accurately during induction and heat 
release than during equilibration. During equilibration all 
other codes tracked the temperature well, with GCKP84 
and EPISODE being superior to the other codes because 
of the use of a much smaller error tolerance. These 
temperature differences show that even for solutions 
generated with large values of EPS, the differences be- 
tween temperature methods A and B were insignificant. 
Furthermore, for test problem 1, solutions generated 
with EPISODE-A and DASCRU-A were as accurate as 
those generated, respectively, with EPISODE-B and 
DASCRU-B, although ERMAX was much larger than 
EPS (table XI). Similar remarks apply to test problem 2 
and the runs with EPISODE and DASCRU. 

The numbers in parentheses in tables XIII and XIV are 
values for the percent differences obtained by ignoring all 
species with standard solution values for mole fractions 
less than 10 -'7. This was done because the selection 
criterion used in the section Efficiency Comparisons 
involved only species with mole fractions greater than or 
equal to 10-7. All mole fractions were larger than 10-7 


for times t greater than or equal to 5 x 10-5 and 10-4 for 
test problems 1 and 2, respectively. 

The mole fraction differences from the standard 
solutions were largest at times immediately after ignition 
and during heat release because the solutions were 
changing rapidly in these regimes. During ignition and 
heat release GCKP84 and EPISODE were inferior to the 
other codes. However, at longer times they were superior 
to the other codes. These results for species mole 
fractions are consistent with the solutions for the 
temperature. In addition, as discussed in the section 
Efficiency Comparisons, although the use of larger error 
tolerances gave poor solutions during ignition and heat 
release, the predicted solutions were satisfactory at longer 
times. Hence, where the main objective of modeling is to 
predict pollutant (e.g., NOx) formation, the user can 
specify fairly large error tolerances without incurring 
severe error penalties. Again, we note that the differences 
between the solutions generated with methods A and B 
(and the same value of EPS) are not significant. Among 
the codes that use a large error tolerance, CHEMEQ was 
the most accurate during ignition and early heat release; 
during late heat release and equilibration, LSODE and 
CREKID were superior to the other codes. 

EPISODE was inferior to LSODE and CREKID at 
early times because of the difference in the error control 
performed by these codes. In contrast to the other codes 
for which EPS is the local relative error tolerance, EPS is 
a mixed error tolerance for EPISODE; that is, it is 
relative for species with nonzero initial mole numbers and 
absolute otherwise. Since most of the species had zero 
initial mole numbers for both problems examined in this 
study (tables IV and V), the error control performed by 
EPISODE was mostly absolute. For pure relative error 
control the estimated local truncation error Ej in species i 
satisfies 

Ei^^v^xrii (17) 

where €rei is the local relative error tolerance. For pure 
absolute error control Ei satisfies 

^/ — ^abs ( 18 ) 

where Cabs is the local absolute error tolerance. Equations 
(17) and (18) show that, since /i, < < 1, relative error 
control is more accurate for the same value of the local 
error tolerance. In other words, to attain comparable 
accuracy for mole numbers that were initially zero, 
EPISODE requires lower values of EPS than the other 
codes. Equations (17) and (18) also show that the error 
controls are of comparable accuracy for tij of order 
0(€abs/€rei)> ^ result that can be used to examine the 
difference in accuracy attained by EPISODE and 
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TABLE XIII.— DIFFERENCES FROM STANDARD SOLUTIONS FOR 
TEST PROBLEM 1 


[Numbers in parentheses are values obtained by ignoring species with mole fractions less than 
10“L] 


Method 

Time, /, s 

9xl0-« 

10-5 

5x10-5 

10-“' 

5x10-'' 

10-3 

Difference in mole fractions, (rms), percent 

GCKP84 

4.001 X 103 

3.866x10' 

2.076 

0.795 

0.138 

0.121 


(266,0) 

(228.8) 





CREKID 

689.5 

2.822x103 

.502 

.515 

.274 

.257 


(3.965) 

(4.008) 





LSODE-A 

766,6 

3.082x103 

.619 

.700 

.512 

.122 


(10.05) 

(9.280) 





LSODE-B 

2.002 xlO^ 

1.044 xlO'O 

.871 

.764 

.335 

.184 


(10.32) 

(9.687) 





EPISODE-A 

4.781 xlO^ 

5.659x104 

2,157 

.784 

.0873 

.0701 


(310.4) 

(267.0) 





EPISODE-B 

7.702x104 

6.789x10’ 

2.182 

.791 

.0881 

.0725 


(317.2) 

(273.0) 





CHEMEQ-A 

647.1 

2.618x103 

10.48 

15.64 

32.72 

42.17 


(0.305) 

(0.299) 





CHEMEQ-B 

649.4 

2.622x103 

5.306 

10.39 

19,17 

20.80 


(0.154) 

(0.150) 





DASCRU-A 

1.785x103 

9.228x103 

1.174 

.435 

.225 

.350 


(107.8) 

(94.66) 





DASCRU-B 

1.819x103 

9.486x103 

1.183 

.430 

.232 

.391 


(110.8) 

(97.14) 





Difference in mole fractions (maximum), percent 

GCKP84 

1.324x10^ 

1.282x105 

5.528 

1.957 

0.451 

0.360 


(374.5) 

(330.7) 





CREKID 

2.285 X 103 

9.359x103 

-.862 

.893 

-.584 

.598 


(5.485) 

(5.746) 





LSODE-A 

2.541X103 

1.022x104 

-1.398 

-1,498 

-1.033 

-.234 


(13.72) 

(13.10) 





LSODE-B 

6.641 X 10^ 

3.464 xlO'O 

-1.695 

-1.624 

-.697 

.459 


(14.42) 

(14.25) 





EPISODE-A 

1.583x10' 

1.877x105 

5.749 

1.895 

.285 

.209 


(439.6) 

(386.5) 





EPISODE-B 

2.555x103 

2.252x10* 

5,814 

1.917 

.288 

.212 


(449.6) 

(395.1) 





CHEMEQ-A 

2.144x103 

8.683 X 103 

19.08 

28.49 

-56.32 

-72.11 


(-0.446) 

(-0.440) 





CHEMEQ-B 

2.152X103 

8.695 X 103 

8.920 

18.32 

36.89 

40.45 


(-0.236) 

(-0.208) 





DASCRU-A 

5.911x103 

3.060x10^ 

3.082 

.953 

-.397 

-.595 


(148.8) 

(135.2) 





DASCRU-B 

6.022X103 

3.146x10^ 

3.114 

.953 

-.424 

-.684 


(153.0) 

(138.8) 





Difference in temperatures, percent 

GCKP84 

0.899 

4.572 

0.208 

0.0795 

0 

0.0381 

CREKID 

.0497 

.111 

.0606 

.0357 

.0333 

.0322 

LSODE-A 

.0587 

.172 

.0951 

.0840 

.0448 

.0360 

LSODE-B 

.0588 

.177 

.136 

.103 

.0501 

.0442 

EPISODE-A 

1.107 

5.548 

.204 

,0628 

.00629 

.0236 

EPISODE-B 

1.137 

5.697 

.207 

.0634 

,00620 

.0236 

CHEMEQ-A 

.0365 

.0430 

-1.466 

-1.635 

-1.917 

-2.163 

CHEMEQ-B 

.0333 

.0409 

.449 

.904 

1.491 

1.539 

DASCRU-A 

.300 

1.590 

.0831 

.0398 

-.0381 

0 

DASCRU-B 

.300 

1.590 

.125 

.0398 

- .0381 

-.0381 




TABLE XIV.— DIFFERENCES FROM STANDARD SOLUTIONS FOR TEST PROBLEM 2 


[Numbers in parentheses are values obtained by ignoring species with mole fractions less than 10"'^.] 


Method 

Time, s 


3xl0-« 

5xl0-« 

10-5 

5x10-* 

10-“' 

5X10-'' 

10-3 

Difference in mole fractions (rms), percent 

GCKP84 

143.6 

8.772x10’ 

8.447 

0.731 

0.615 

0.852 

0.835 


(103.2) 

(15.98) 

(5.925) 

(0.723) 




CREKID 

4.106 

3.570x10’ 

0.322 

0.0458 

.0459 

.213 

.297 


(1.089) 

(0.504) 

(0.222) 

(0.0435) 




LSODE-A 

11.62 

4.486x10’ 

10.62 

0.444 

1.011 

.0721 

.170 


(8.642) 

(1.977) 

(8.086) 

(0.380) 




LSODE-B 

653.9 

10.85 

12.45 

1.854 

2.159 

.563 

.592 


(30.83) 

(7.167) 

(8.741) 

(1.584) 




EPISODE-A 

105.7 

7.612x10’ 

6.996 

0.798 

.276 

.0358 

.0449 


(76.26) 

(14.11) 

(4.806) 

(0.621) 




EPISODE-B 

157.3 

8.743x10’ 

8.542 

0.947 

.358 

.0260 

.0174 


(101.0) 

(15.87) 

(5.883) 

(0.737) 




CHEMEQ-A 

3.830 

3.457x10’ 

0.0528 

2.885 

10.25 

26.71 

37.85 


(0.198) 

(0.0893) 

(0.0382) 

(2.570) 




CHEMEQ-B 

3.839 

3.447x10’ 

0.103 

1.591 

5.990 

25.40 

38.50 


(0.234) 

(0.138) 

(0.0735) 

(1.433) 




DASCRU-A 

30.17 

10.23 

25.90 

0.362 

8.571 

All 

9.640 


(23.24) 

(7.541) 

(1.829) 

(0.280) 




DASCRU-B 

30.18 

1.247x10’ 

2.615 

0.261 

.564 

.641 

7.164 


(23.24) 

(7.085) 

(1.787) 

(0.201) 






Difference in 

mole fractions (maximum) 

percent 



GCKP84 

317.5 

3.397x10* 

21.03 

1.471 

-1.396 

-1.622 

-1.764 


(164.2) 

(-26.75) 

(21.03) 

(1.471) 




CREKID 

- 10.382 

1.383x10* 

0.768 

-0.101 

-.122 

-.446 

-.609 


(1.728) 

(-0.854) 

(0.768) 

(-0.101) 




LSODE-A 

21.85 

1.738x10* 

29.12 

0.965 

-1.706 

.144 

.382 


(13.64) 

(-3.207) 

(29.12) 

(0.934) 




LSODE-B 

1.987x103 

30.55 

36.45 

4.051 

3.932 

1.229 

1.203 


(48.95) 

(-12.25) 

(31.51) 

(3.821) 




EPISODE-A 

222.4 

2.948x10* 

17.13 

2.034 

.717 

.0684 

.131 


(121.5) 

,(-24.06) 

(17.13) 

(1.820) 




EPISODE-B 

309.8 

'3.386x10* 

20.97 

2.422 

.898 

.0621 

.0303 


(160.7) 

(-26.87) 

(20.97) 

(2.179) 




CHEMEQ-A 

- 12.32 

1.339x10* 

-0.151 

-6.203 

-20.52 

-46.38 

-63.40 


(-0.334) 

(0.177) 

(-0.129) 

(-6.203) 




CHEMEQ-B 

-12.51 

1.335x10* 

-0.281 

-3.034 

11.18 

44.45 

65.14 


(-0.402) 

(0.256) 

(-0.260) 

(-2.731) 




DASCRU-A 

57.47 

23.44 

- 100.0 

0.934 

33.19 

.957 

37.25 


(38.03) 

(-14.62) 

(6.530) 

(0.840) 




DASCRU-B 

57.50 

4.830X10’ 

6.359 

0.674 

-2.174 

-1.244 

27.24 


(38.03) 

(-11.89) 

(6.359) 

(0.630) 




Difference in temperatures, percent 

GCKP84 

0.333 

2.203 

0.733 

0.0795 

0 

0 

0 

CREKID 

.00477 

.0875 

.0152 

-.00135 

.0218 

.0126 

-.00641 

LSODE-A 

.0306 

.370 

-.266 

.146 

.203 

- .0179 

- .0282 

LSODE-B 

.105 

1.017 

-.140 

.310 

.373 

.0704 

.0406 

EPISODE-A 

.261 

1.856 

.635 

.0739 

.0361 

- .00448 

- .00292 

EPISODE-B 

.344 

2.194 

.770 

.0910 

.0466 

- .00447 

-.00421 

CHEMEQ-A 

.00486 

.0209 

- .0227 

-.524 

-1.381 

-1.025 

-1.227 

CHEMEQ-B 

-.00189 

.0108 

- .0288 

-.107 

.0265 

1.283 

1.678 

DASCRU-A 

.0665 

.796 

.262 

.0398 

.0367 

0 

.0687 

DASCRU-B 

.0665 

.796 

.262 

0 

0 

- .0687 

-.138 




LSODE, for mole numbers that are initially zero. Using 
the values for the error tolerances given in tables XI and 
XII, we found that the two codes were of comparable 
accuracy for «,• of order 0(10-4) and 0(10-3) for test 
problems 1 and 2, respectively. It should be emphasized 
that these values are only estimates because the two codes 
control the root-mean-square norm of the estimated local 
errors and not the estimated local error for each species. 

These estimates for mole numbers correspond to mole 
fractions of order 0(10-3) and 0(10-2) for test 
problems 1 and 2, respectively, assuming a mixture mean 
molar mass of order 0(10). Tables IV and V and figures 1 
and 2 show that these values were not attained at early 
times by species with initial zero mole numbers. Hence 
EPISODE is expected to be inferior to LSODE during 
ignition and early heat release. However, during late heat 
release and equilibration, when many of the mole 
numbers exceed these values, EPISODE is expected to be 
superior to LSODE. Examination of tables XIII and XIV 


confirms this behavior. We note that, relative to LSODE, 
EPISODE was (1) inferior at short times, (2) superior at 
long times, and (3) comparable at intermediate times. 

These observations imply that the error control used by 
EPISODE is unsatisfactory for problems of the type 
examined in this study. Many of the variables have zero 
initial values and they never reach unity. The nature of 
the error control performed by EPISODE therefore 
requires small values of EPS for acceptable accuracy at 
short times. The continued use of these small values of 
EPS at long times is wasteful because it results in 
solutions that are more accurate than required by the 
selection criterion. A simple method for increasing the 
efficiency of EPISODE is therefore to switch to pure 
relative error control (lERROR = 2) with a larger value of 
EPS once the species mole numbers have reached 
acceptable values. Tables XV and XVI present the effects 
of such a switch. In these tables, /switch is the time at 
which the switch was made. The program was run up to 


TABLE XV.— EFFECT OF SWITCH TO PURE RELATIVE ERROR 
CONTROL ON COMPUTATIONAL WORK REQUIRED FOR 
TEST PROBLEM 1 


(a) EPISODE-A 


Time at 
which 
error 
control is 
switched, 

^switch » 

S 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number 

of 

functional 

evaluations, 

NFE 

Total 
number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total 

number of 
reaction 
rate 

constant 

evaluations, 

NRRC 

Total 

CPU 

time 

required, 

CPU, 

s 

10-5 

121 

223 


175 

0.37 

1.5x10-5 

150 

290 


273 

oo 

2x10-5 

164 

314 

40 

294 

.51 

2.5x10-5 

175 

332 

40 

304 

.53 

5x10-5 

194 

369 

42 

347 

.58 

lo--* 

209 

394 

41 

367 

.62 

5X10-'* 

236 

444 

44 

415 

.68 

aiO-3 

244 

463 

41 

435 

.70 


(b) EPISODE-B 


Time at 
which 
error 
control is 
switched, 

^switch* 

s 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number 

of 

functional 

evaluations, 

NFE 

Total 
number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total 

number of 
reaction 
rate 

constant 

evaluations, 

NRRC 

Total 

CPU 

time 

required, 

CPU, 

s 

1.5x10-5 

140 

257 

37 

240 

0.43 

N> 

X 

O 

1 

148 

278 

37 

258 

.45 

2.5x10-5 

157 

292 

36 

277 

,47 

5x10-5 

181 

340 

37 

322 

.53 

lo--* 

201 

367 

39 

350 

.57 

o 

X 

234 

436 

41 

413 

.66 

aio-3 

248 

460 

36 

447 

.68 


^No switching performed. 


19 









TABLE XVI.— EFFECT OF SWITCH TO PURE RELATIVE ERROR 
CONTROL ON COMPUTATIONAL WORK REQUIRED FOR 
TEST PROBLEM 2 


(a) EPISODE-A 


Time at 
which 
error 
control is 
switched, 

^switch* 

S 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number 

of 

functional 

evaluations, 

NFE 

Total 
number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total 

number of 
reaction 
rate 

constant 

evaluations, 

NRRC 

Total 

CPU 

time 

required, 

CPU, 

s 

3x10-6 

92 

167 

35 

132 

0.53 

4x10-6 

119 

212 

35 

172 

.64 

5x10-6 

101 

186 

31 

168 

.58 

10-5 

112 

204 

36 

188 

.65 

X 

o 

( 

118 

218 

34 

204 

.68 

10-4 

125 

227 

35 

215 

.72 

5x10-6 

128 

233 

35 

219 

.72 

aiO-5 

127 

227 

31 

229 

.71 


(b) EPISODE-B 


Time at 
which 
error 
control is 
switched, 

^switch* 

S 

Total 
number 
of steps 
required, 
NSTEP 

Total 

number 

of 

functional 

evaluations, 

NFE 

Total 
number 
of Jacobian 
matrix 
evaluations, 
NJE 

Total 

number of 
reaction 
rate 

constant 

evaluations, 

NRRC 

Total 

CPU 

time 

required, 

CPU, 

s 

4x10-6 

88 

173 

30 

140 

0.54 

5x10-6 

104 

215 

32 

170 

.62 

10-5 

119 

221 

36 

191 

.67 

Uh 

X 

O 

<-* 

127 

242 

32 

219 

.71 

10-6 

130 

248 

33 

229 

.73 

5x 10-4 

142 

296 

37 

258 

.84 

aiO-5 

145 

303 

34 

273 

.86 


^No switching performed. 


time t = tswitch with IERROR = 3 and EPS of 10-6 and 
10-5 for test problems 1 and 2, respectively. At t = ts witch 
the integrator was reinitialized, the error control was 
switched to IERROR = 2, and EPS was increased to 
10-2. The problem was then run to completion with these 
new values for lERROR and EPS. Note the significant 
decreases in the computational work obtained by 
switching to pure relative error control. The resultant 
CPU times compare favorably with those required by 
LSODE. This switching process is, however, unsat- 
isfactory for the following reasons: 

(1) The value of ^switch that minimizes the CPU time is 
a function of the problem and the temperature method 
used. Moreover, a poor choice for ^switch can make the 
run prohibitively expensive. For example, for test 
problem 1 the run with EPISODE-B and fswitch = 10“^ 
was not successfully completed even after a CPU time of 
2 min. 

(2) The switching requires unnecessary computational 


work associated with the reinitialization and restart of the 
integrator. 

(3) Using low values of EPS for /switch results in 
tighter error control than is necessary to satisfy the 
accuracy criterion for species with initially nonzero mole 
numbers and for the temperature. 

A more appropriate error control for problems of the 
type examined in this study is therefore one that is 
absolute for mole numbers that are initially zero and that 
automatically becomes relative once these mole numbers 
reach acceptable values. The error control should also be 
relative for species with initially nonzero mole numbers 
and for the temperature. Such an error criterion can be 
realized by requiring Ej to satisfy 

Ei S ereltt/ + Cabs,/ (19) 

where Cabs,/ is the local absolute error tolerance for 
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Time, t, s 


(a)GCKP84. (b) CREKID. (c) EPISODE. (d) LSODE. 

Figure 5. — Variation with time of the step length successfully used by GCKP84, CREKID, EPISODE, and LSODE for test problem 1. 


species i. For values of «,■ < (eabs,/%el) the error control 
is mostly absolute, and for rij > (eabs./crel) it is mostly 
relative. This is, in fact, the nature of the error control 
used in LSODE; it is therefore more efficient than 
EPISODE for integrating combustion kinetic rate 
equations. 

Step Length Comparisons 

All codes used in the present study automatically select 
a step length during the course of the integration. Some 
of the codes (GCKP84, DASCRU, and EPISODE) 
require trying a user -supplied initial value. The other 
codes automatically select the value for the initial step 
length. The size of the step successfully used by the code 
indicates both the efficiency of the code and regions 
where difficulties due to stiffness arise. Figures 5 to 8 
present plots of the step length used by each code through 
the course of each problem. To facilitate comparisons 
(among the faster codes) at early times for test problem 2, 
figure 9 presents the variation of the step length between 
t = 10-6 and t= 10-5. 


Figures 6 and 8 illustrate the small steps that classical 
methods have to use to ensure solution stability. For both 
test problems the explicit Runge-Kutta technique used 
small step lengths to track the solutions through 
induction and heat release. During equilibration the step 
lengths remained small, and thus prohibitive amounts of 
computer time were required. The difficulties with 
CHEMEQ (figs. 6 and 8) included the selection of a very 
small initial step length, the continued use of small step 
lengths because of the very small increases allowed after 
satisfactory convergence, and its inability to select a 
suitable step length during equilibration. Much computer 
time was wasted in the search for an appropriate step 
length. In addition, the search was restricted to very small 
step lengths. These factors make CHEMEQ very 
expensive to use. 

We note that all codes use small steps during induction 
and early heat release. In these regimes the species and 
temperature change rapidly (figs. 1 and 2). Most of the 
species and temperature have positive time constants — an 
indication that the differential equations are unstable. 
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Hence the step lengths are constrained to small values. 

For test problem 1 , EPISODE was superior to 
GCKP84 and CREKID was superior to LSODE. Note 
that although in the postinduction regime EPISODE 
selected steps comparable to those selected by LSODE, 
its difficulty in tracking the solution during ignition and 
early heat release made it less efficient than LSODE. For 
test problem 2, at long times, CREKID took step lengths 
comparable to those of LSODE. However, at times 
preceding and immediately after ignition (t^3 /^s), 
CREKID took much smaller step lengths and was hence 
less efficient than LSODE. Although not presented 
herein, the run with CREKID and EPS = 10-2 showed 
that this code had difficulty selecting a suitable step 
length at early times. Much computer time was wasted by 
repeatedly unsuccessful attempts at selecting a larger step 
length. This run was therefore less efficient than the one 
with EPS = 10-3 (fig. 4). The step lengths selected by 
EPISODE were comparable to those selected by 
GCKP84, although the former used larger steps initially. 

The results discussed above indicate that the step 
length to be used is regime dependent: during induction 
and early heat release, when the solution changes rapidly, 
small steps have to be taken to ensure solution stability. 
The search for large step lengths is both futile and 
expensive. At later times, however, when the differential 
equations are more stable, larger step lengths can be 
used. In these regimes it is worth attempting to use a 
larger step length after every step. 


Summary of Results and Conclusions 

Several algorithms (GCKP84, CREKID, LSODE, 
EPISODE, CHEMEQ, and DASCRU) for numerically 
integrating stiff ordinary differential equations arising in 
combustion chemistry have been compared. The 
following performance indicators were recorded for 
purposes of this comparison: 

(1) Total number of steps required to solve a problem 

(2) Total number of derivative evaluations 

(3) Total number of Jacobian matrix evaluations 

(4) Total number of rate constant evaluations 

(5) Total CPU time required to solve a problem 

In addition, the errors relative to reference solutions 
(generated with LSODE with a very low relative error 
tolerance) were recorded at selected times. These tests 
were conducted on two combustion kinetics problems: 
one involving 11 species and temperature with 12 reac- 
tions, and the other involving 15 species and temperature 
with 30 reactions. Both problems included all three 
combustion regimes: induction, heat release, and 
equilibration. 


This study has shown that for comparable accuracy the 
fastest package currently available for integrating 
combustion kinetic rate equations is LSODE. This merits 
special notice because LSODE was developed as a 
multipurpose stiff differential equation solver with no 
one particular application as its objective. Its disadvan- 
tages, however, include a large storage requirement and a 
relatively long startup time. Where storage is a significant 
problem, both EPISODE and CREKID are attractive 
alternatives. However, a poor guess for the initial step 
length to be tried by the integrator can make EPISODE 
prohibitively expensive to use. It can also result in 
incorrect and unstable solutions. Some experimentation 
with different values for the initial step length may be 
necessary to obtain its optimum value. In addition, the 
error control performed by EPISODE was found to be 
unsatisfactory for problems of the type examined in this 
study. For acceptable accuracy at short times small values 
of the error tolerance had to be used. Significant 
decreases in computational work were obtained by 
switching to pure relative error control in the 
postinduction regimes. An error control that is more 
appropriate for combustion kinetics problems is required 
for additional increases in computational efficiency. The 
code CREKID needs refinement in the area of step- 
length selection, especially at early times, before 
significant increases in its speed can be realized. 

The step length to be used in integrating combustion 
kinetic rate equations was found to be regime dependent. 
During induction and early heat release, small steps had 
to be taken to ensure solution stability. At later times, 
however, when the differential equations were more 
stable, larger steps could be used. 

An important conclusion from this study is that using 
an algebraic energy conservation equation to calculate 
the temperature does not result in significant errors or 
inefficiencies. On the contrary, this method can be more 
accurate and efficient than evaluating the temperature by 
integrating its time derivative. 

Another important conclusion is that, where the 
objective of modeling is to predict pollutant (e.g., NOx) 
formation, the user can specify large error tolerances 
without incurring significant error penalties. 

A simple method for realizing significant efficiency 
increases was demonstrated. This involved updating the 
reaction rate constants only for temperature changes 
greater than an amount that was problem dependent. An 
expression for the maximum temperature change allowed 
before such an update was presented and shown to result 
in significant savings. 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, June 6, 1984 
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Appendix A 

Outline of Methods Studied 


The ordinary differential equations (2) and (8) 
describing homogeneous gas-phase chemical reactions 
can be generalized as follows: 


yi=^‘=fi(yk) i,k=l,N 
yi(t = 0)=yi(0)=&ven 


(Al) 


where for temperature method A (see the section 
Evaluation of Temperature) 


yi=rii / = 1,NS 

7V=NS 


and for temperature method B 


(A2) 


this study, some (CHEMEQ, CREKID, and DASCRU) 
are single-step methods; the rest (EPISODE, LSODE, 
and GCKP84) are multistep methods. A single-step 
method provides a rule for computing y„ at t„ from a 
knowledge of the step length and the solution i at 
the previous step t„-\. A multistep method, on the other 
hand, uses and the solution at several earlier points to 
generate y„. For example, the ^-step method uses h„ and 
the K earlier values Jn - 1 , - 2 . • • • , - k (at times 

tn-2» •••» tn-id to compute Since at the outset only 
the initial conditions XO) are known, multistep methods 
are not self-starting. EPISODE and LSODE resolve this 
difficulty at the initial point by starting with a single-step, 
first-order method. As the integration proceeds, the 
solutions that they generate provide the necessary values 
for a multistep method. 

We now outline the methods studied in the present 
work. In particular, we examine how each method 
advances the solution by one step (of length h„) from 
time _ 1 to time t„. To avoid confusion between the time 
step length h„ and the molal-specific enthalpy for species 
«, the latter is denoted by fi„. 


yi=rij /=1,NS 
.VnS +1 = T 
7V=NS+1 


In vector notation, equation (Al) becomes 


EPISODE and LSODE 

Both EPISODE and LSODE use the integration 
(A3) formulas developed by Gear (refs. 20 and 21). These 
formulas involve linear multistep methods of the form 


y=i j=0 


r (A4) 

^(/ = 0)=j(0)= given ) 


where the underscore is used to denote a vector quantity. 
A matrix is denoted by a boldface letter. This notation is 
used throughout this appendix. 

The initial-value problem is to determine values for yj 
(/= 1,A/) at the end of a prescribed time interval, given // 
(/= l,N) and the initial values /,• (0) (/ = l,N). 

All methods examined in the present study are step-by- 
step methods. They compute approximations y„{=yi^„; 
i=l,N) to the exact solution X/„) at discrete points tj, t 2 , 
ti, ... The size h„ of the spacing ( = r„-/„_i) may vary 
from one step to the next. Of the methods examined in 


where yj„ is an approximation to the exact solution 
yt (fn)y yi,n\. =fi is an approximation to the exact 
derivative y, (t„){=// b>k (/„)]], and the a„j and 

are associated with the particular formula 
selected by the user. The options include implicit Adams 
predictor-corrector formulas and backward differen- 
tiation formulas (BDF). As discussed in appendix B, 
BDF’s were found to be more efficient for problems 
examined in this study. The discussion is therefore 
restricted to BDF’s. For a BDF of order q, K\ = q,K 2 = 0, 
and equation (A5) reduces to 

(1 

yi,n~ 2 kfifffi^Qyj^fi (A6) 

y=i 

The step length = can vary from one step 

to the next in EPISODE but is held constant for q+\ 
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consecutive successful steps in LSODE. Hence, for 
EPISODE, a„j and can vary from one step to the 
next, but in LSODE they are predetermined constants 
corresponding to the order used. The predicted values of 
y„ and y„, denoted by and respectively, satisfy 
the equation 


yTn~'^rt^n.oyfl= ««,>)'/.« -7 


so that equation (A6) can be written as 


=yf„ - ^ A.ai'S * A.Qy».n 


(A7) 


At each integration step equation (A7) must be solved 
for y„. This is accomplished as follows: let 

G(y„)=yft — hff0n,(^n~ (A8) 


For this method, much computation time is involved in 
forming the Jacobian matrix J and in doing the linear 
algebra necessary to solve equation (A9). To reduce this 
computational work, the matrix P is not updated at every 
iteration. For further savings, it is updated only when it 
has been determined to be absolutely necessary for 
convergence. Hence the iteration matrix is only accurate 
enough for the iteration to converge, and the codes may 
use the same matrix over several steps of the integration. 
In any case, both EPISODE and LSODE update P at 
least every twentieth step. 

The predicted values yfl and yfl are required before 
equation (A9) can be solved for These are obtained by 
a Taylor series expansion as follows: the history of the 
solution is maintained in the Nordsieck array (which is a 
Taylor series array) z„ of size Nx (g + 1). The /th row Zi,n 
contains the g + 1 elements 

- • Hfi •• frti 

Zi,n~yi,n> ^nyi,n> •••> 


then the solution to equation (A7) is equivalent to finding 
the zero of G. Equation (A8) can be solved in a variety 
of ways. For stiff problems the most efficient is a variant 
of Newton’s method: 


where y^], is the approximation to iefyj/d^^. Thus, if 
z„-i has been obtained. 



(A9) 

where m+\ is the current iteration number and the 
matrix P is given by 


where is given by 



j<k 


] 


y. At — 0,1,2, 




(AlO) 


=«ik-M«.o/iiS 

where djic is the Kronecker symbol, 

. r= 0 iVA: 

I (=* 

a«d 4".’ = l,N) are elements of the Jacobian 

matrix j’ 


where is the binomial coefEcient, defined as 
(^) " klU-k)l 

Thus the predicted values are obtained by a simple 
Taylor series expansion by using equation (AlO). Once 
the predicted values have been calculated, equation (A9) 
is iterated until convergence. The test for convergence of 
the iterates y^^^is based on successive differences of these 
quantities and the local error tolerance EPS as discussed 
below. For EPISODE a vector YMAX is constructed as 
follows: 



IERROR= 1 (absolute error control): 
YMAX,= 1 1 = 1,2 N 
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IERROR=2 (pure relative error control): 


EPISODE: 


lERROR = 3 (semirelative error control): 

YMAXj = max ([>'/,„ - 1 1, - 2 1) for i that satisfy yi(0) ^ 0 

= max {1 , \yj^„ _ 1 1) for / that satisfy ^,(0) = 0 

Convergence is said to occur if 

^e'=C£EPS (All) 

where e' is proportional to EPS, the local error tolerance. 

For LSODE and option ITOL = 2 (see appendix B for 
other options included in this package), an error vector 
EWT is constructed as follows: 



EWT,-— €re]|y^n_l| + eabs,/ /— 1,2,...,A/^ 

where erei(=EPS) and eabs,/ ar® the local relative error 
tolerance and the local absolute error tolerance, 
respectively, for the variable /. 

Convergence is said to occur if 


sci (A12) 

In equations (All) and (A12) the proportionality 
constants cg and Cg are chosen to make the convergence 
tests consistent with the local truncation error. If 
convergence does not occur after the Erst iteration, both 
EPISODE and LSODE anticipate the magnitude of dm 
one iteration in advance by assuming that the iterates 
converge linearly. Thus dm+\> which does not yet exist, is 
estimated by 



s/£EPS 


LSODE: 




where tg and tg are test constants that depend on the 
order used. If the error test fails, the step is repeated with 
a different h„ or a lower order until either the test is 
passed or the situation is considered hopeless. If the test 
is passed, the step is accepted as successful, and the entire 
Nordsieck array is updated by 

7 — " 4 - ^ / 

*•11 *‘n ^EnLn 


where 




and the row vector ^„(/i^„;/=0,l g) is determined by 

the formulas used and satisfies /o,« = 1 and /i,/i = l//3„,o* 
For EPISODE, /„ depends on the variable step length 
and is computed at the start of every step. For LSODE, 
are constants that are set at the beginning of the 
problem. 

In summary, the predictor and corrector steps are 
given by 

Predictor: 





Corrector: 


If the corrector iteration fails to converge in three 
iterations, h„ is reduced if P is current and the step is 
retried; otherwise, P is updated and the step is retried. If 
the corrector converges after M (Ms 3) iterations, an 
estimate of the truncation error is made and is accepted if 
it passes the following test: 


= -G\y„('»)] > /n = 0,l,...,M-l 
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and 


M-\ 

£«= E 

w = 0 

z„=z<g> + e^ 


NS 

E Qk/r* 

Dj-= — 

NS 

yk^p,k 


The objective of this decomposition is to enable 
CHEMEQ factorization of from Dj 

In this technique, developed by Boris and Young (ref. 

11), the rate equation (A2) is expressed as a difference 

between two positive-definite terms as follows: ^i=^iyi= - 


^£=fi=Qi-Di (A13) 

where, for species i, the production rate Qi and the 
destruction rate can be derived from equation (3) 

JJ \ 

Qi=p~^ E J 

j=i I 

> (AM) 

JJ I 

E (piiRj+njR-j) ) 

j=i ' 


When temperature is an independent variable (method 
B), its rate equation (eq. (8)) can be cast in a similar form 
by combining equations (8) and (A 13) 


where Lj is obtained simply by dividing £),• by y (i.e., 
Li=D/yj). With this new notation, equation (A13) can 
be written as 

^'=G, = (2,-^ (A15) 

which, for constant Q, and L,-, can be solved to give 

yii^n) -yiifn - 1 + = ^ + ^iifn “ ~ ^ J 

xexp( -£,/»„) (A16) 


dt dt 


NS 

“ E fkFk 

k=l 

NS 

^Ej ykpp.k 


NS 

- E (^Qk-Dk)fTk 

k=\ 

NS 

E yk^p.k 

k=\ 


With this expression, it can be seen that 1/L,( = t/) 
describes how quickly the variable y( reaches its 
equilibrium value. 

In advancing the solution from time /«_ i to time t„, all 
of the equations are separated into two classes, stiff and 
nonstiff, according to the criterion 


= Qns+i-T>ns+i 


= — Dp 


h„ r>l stiff 

T/„_j |^<1 nonstiff 


where 


Qt= 


NS 

E 

k=\ 

NS 

E y/^p.k 

k=V 


where t/ „ _ j denotes the value of t,- at time - 1 • The two 
types of equations are integrated by separate predictor- 
corrector schemes. For equations classified as nonstiff, 
the improved Euler method— the Euler method as 
predictor and the modified Euler method (or trapezoidal 
rule) as corrector — is used. For equations classified as 
stiff a simple asymptotic formula is used. 
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Predictor (nonstiff): 


Hnfi.ti - 1 


Predictor (stiff): ( (A 17) 


y(0) _ yj.n - l(^'^t.n - 1 ^n) + - lQi,n - 1 


Corrector (nonstiff): 

Corrector (stiff): \(A18) 

[y [*>+'«-.] [e!5> +Gu-i] 

- 1 - 1 ^nj ^ ^ 

In equations (A17) and (A18), m + 1 is the current 
iteration number. The zeroth iterate is the result of the 
predictor step. Also, Convergence is 

ascertained by comparing jy ^ ^ with for all N 
equations using the relative error criterion 


vJw + l) _v(w) 




<EPS 


(A19) 


min 


If after ITMAX iterations any of the variables fails 
the test given by equation (A19), the step length is halved 
and the step repeated. If all AT variables pass the test after 
M iterations (M< ITMAX), the step is accepted as 
successful and the solution is updated 


code consists of two algorithms developed for the two 
distinctly different regimes identified in the section 
Accuracy Comparisons. These regimes are (a) induction 
and early heat release, when the ODE’s are dominated by 
positive time constants, and (b) late heat release and 
equilibration, when the ODE’s are more stable. Both 
algorithms are based on an exponentially fitted 
trapezoidal rule, but they use different iterative methods 
for convergence. 

In the CREKID method the temperature is not treated 
as an additional independent variable, so the number of 
ODE’s is equal to NS. The temperature is calculated from 
the algebraic energy conservation equation (eq. (7)). In 
the following discussion the variables y,- (i = l,NS) 
therefore refer only to the species mole numbers. 

The species rate expression/) (the right side of eq. (Al)) 
can be expanded in a first-order-truncated Taylor’s series 
about the current approximate solution, y„_i 

* = 1. NS), as follows: 


/) ~fi,n - 1 + ^ yi,n-l) 

\^yt/ yn-\ 


(A20) 


where dfi/dyi is the total derivative of // with respect to yi 
and is given by the chain rule as 


g a/) dy^dt ^ dfi dT/dt 


dyi 1 ^yk dyi/dt dt dyi/dt 


_ 1 / ^ 

fi\ki 9Tdt) 


To avoid confusion with the partial derivative of /) with 
respect to (dfi/dy^, a special notation is employed 


7^<Mi 

' dyi 

With this new notation equation (A20) becomes 


V. /=1 N 

yi.n yi,n ' 


fi~fi,n-l'^^i,n-l(yi 1) 


(A21) 


or 


CREKID 

In CREKID, attention is paid to the distinguishing 
physical and computational characteristics of the 
induction, heat release, and equilibration regimes. This 


~fi,n - 1 ^i,n - iCV/ yt,n - 1) 
which can be integrated to give 
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yi,n '^y^n - 1 ^ ^nf i,n ~ 1 


~ QXp(hriZi^ri-\)- 




(A22) 


Consider now, an approximate solution to equation (2) 
based on a variation of the second-order trapezoidal rule, 
the tunable trapezoid. 


yi,n~yi,n — \ — ll 


(A23) 


The substitution of equation (A22) into equation (A23) 
gives 




• + 


h„Zi.n-\ l-exp(/i„Z/,„_i) 


(A24) 


molal-specific enthalpy of species k at temperature 
and Hq (Tq) is the initial mixture mass-specific enthalpy 
at the initial temperature Tq. 

Newton-Raphson corrector equations with log variable 
corrections (for self-scaling of the widely varying mole 
numbers) are given by 


NS 
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Alogyg 
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Alog 7’^'”) = 


aiogr("’) 


t,n 


NS 


dF 


(m) 


/t = i dlogy^'^) 




+ 




aiog7('”> 


(A27) 




In order to maintain absolute A-stability of equation 
(A23) (i.e., y^n remains bounded as h„ is increased 
indefinitely), f/, must be restricted to the interval (0.5, 
1.0). For values of Z,->0, equation (A24) gives f/,<0.5. 
CREKID resolves this problem by setting Z/=0 whenever 
it is greater than zero. This gives C/, = 0.5, so that 
equation (A23) defaults to the second-order-accurate 
trapezoidal rule. However, for Z,-s0, equations (A23) 
and (A24) together are equivalent to the locally exact or 
exponential solution, which has an equivalent polynomial 
accuracy of order six to eight. Thus equations (A23) and 
(A24), with the constraint (0.5 < t/,< 1), constitute an 
exponentially fitted trapezoidal rule, a method which is 
A-stable and has a polynomial-order accuracy of at least 
two and as great as six to eight. 

At each integration step, equation (A23) must be 
solved for This is accomplished by Newton-Raphson 
iteration in regime b and Jacobi-Newton iteration in 
regime a. 

A Newton-Raphson functional = l,NS) for the 

species mole numbers is defined from equation (A23) by 


W.n 


i,n 



For temperature the functional is defined from the 
enthalpy conservation equation (7) as 


To start the iteration process, the predicted values yf'^ 
and are obtained quite simply by setting them equal 
to the values at the previous time step 


yfn=yi,n-l /=1.NS 


The Jacobi-Newton iteration technique can be derived 
from the Newton-Raphson iteration procedure by 
neglecting the off-diagonal elements of the Jacobian 
matrix dfj/dy/(. With this simplification, equations (A27) 
reduce to 


Alogy^"") = 

aiogy(/ 


(A28) 

(A29) 


The iteration procedure is further simplified as follows: 
the expression for log , derived from equa- 

tion (A25), is given by 


NS 

45 = E ySfr^m-ffoiTo) (A 26 ) 

k= 1 

where m is the iteration number, is the mth- 

approximation to the exact value r(/„), is the 


= ASL 

dlogyj'H^ hn^i.n dy\^^ 

In deriving an expression for 9f\^^/dy\^\ the partial 
derivatives with respect to the inverse mean molar mass 
are assumed to be negligible in comparison with the 
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other terms. This results in the following approximation 
for (see eq. (11) for the exact expression for 






For Jacobi -Newton iteration, this expression is further 
approximated by 


7’W = 7'„. 


1 + 


NS 

k=l 


NS 

^ y^^rfp.k'Rn-X 
k=\ 


For both Newton-Raphson (NR) and Jacobi-Newton 
(JN) iteration schemes the current values 4^^*^ 
j’(w + 1) updated by the approximate equations 




(yilRj + v(^R-j) 


which when combined with equation (A14) gives 


y\T'^ =4n^[i+Aiog4^>] ) 

> (A30) 

7<"*‘l=7'<"'>[l+Alog7;<”)] \ 




The test for convergence of the iterates is based 

on the value A log and is given by 


where is the destruction rate of species i. With these 
simplifications equation (A28) can be solved explicitly for 
the iterative corrections 


Alog^^^ = - 



From equation (A26) the following expression can be 
derived; 






(A31) 


This test is used only with variables whose magnitudes 
are greater than 10-20; that is, the summation does not 
include species with mole numbers less than or equal to 
10-20. If convergence is not obtained after ITMAX 
iterations, the step length is reduced and the step retried. 
If convergence is achieved in M iterations (Ms ITMAX), 
the step is accepted as successful and the solution is 
updated 


where is the constant pressure molal-specific 

heat of species k at temperature 7’^'”^ . Substituting the 
preceding equation into equation (A29) gives 

~P^n 

Alog = ’1 



To start the iteration process, the predicted values 
are obtained from equation (A22) 


yi,n-\y i,n-\ 


txp{hfjZj fj - ]) 

hnZj f, _ ] 
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The predicted temperature T^^'> is obtained by a single 
Newton-Raphson iteration 


y. =y(A^ 7 = 1 NS 

■'i.n ■'i.n ‘ 
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CREKID automatically selects the iteration scheme 
(JN or NR) to be used for solving equation (A23). During 
induction and heat release, when small step lengths are 
required for solution stability, the JN iteration is used to 
minimize computational work. During late heat release 
and equilibration, when the ODE’s are more stable and 
larger step lengths can be used, NR iteration is preferred 
since it has a much larger radius of convergence than JN 
iteration. The regime identification test exploits the fact 
that during equilibration many reactions achieve a 
condition in which the forward and reverse rates are large 
but have vanishingly small differences. The actual test 
employed at the beginning of each time step is 
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t/;|<10-3(Q^- + Z?/) (A32) 

where Q/ and Z)/ are the production and destruction 
terms, respectively, for species i (eq. (A14)). If any of the 
species satisfies equation (A32), regime b is obtained and 
the NR iteration is used for the step. If none of the 
species satisfies equation (A32), regime a is obtained and 
the JN iteration is used for the step. 

For this method much computation time can be 
involved in calculating Z/ because it requires the 
evaluation of the Jacobian matrix. To avoid evaluating 
this matrix, is estimated as follows: if equation 

(A21) is applied to the entire time step 

fi,n^fi,n-~ \ “A/? - l) 

which, when substituted into equation (A22), gives 

fi,n=hn~\^^V(hnZi^n.{) 


conditions may arise, for example, in multidimensional 
modeling because of the averaging of mole numbers over 
adjacent grid nodes. CREKID therefore ‘‘filters'" the 
initial conditions to provide physically meaningful initial 
mole numbers and net species production rates. For 
purposes of this filtering CREKID uses the 
decomposition performed in CHEMEQ (eqs. (A 13) and 
(A 15)). On the first call to CREKID it uses this 
formulation over one time step of length h\, which is 
determined by 


/ 

The predictor-corrector algorithm uses equation (A 15) 
as the predictor 

An implicit Euler corrector is then iterated to 
convergence 




T-lOg 

fin 



However, fj „ is not known at the start of the step, so 
approximations have to be developed for Z,' j. For the 
NR iteration, Z/_„_i is approximated by 


T-log 

rin 



for <1 

fi,n-2 


- 0 otherwise 


In these equations 7,(0) are the initial values and the 
subscript 1 is used to indicate that this is the first step. 
Using equations (AI3), (AI5), and (A30), together 
with the approximations and + 1) = 

^ , the preceding corrector equation can be rewritten 
to provide the following expression for the log variable 
corrections A log 


Alog^^f ^ = 


yif + hiDif 


(A33) 


For the JN iteration, Z,- „„i is approximated by 


Z _ ^ fi,n- 

(\n-l- T 


fi.n - 
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Equation (A33) is iterated until convergence, which is 
ascertained by the criterion given by equation (A31). If 
convergence is not obtained after 10 iterations, the step 
length is halved and the step retried. If convergence is 
obtained after A/ iterations (M<10), the step is accepted 
as successful, the solution for the mole numbers is 
updated 


= 0 otherwise 


CREKID also includes an algorithm for filtering the 
initial conditions that may be ill posed. These ill-posed 


and the temperature T\ is obtained by a single Newton- 
Raphson iteration 
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GCKP84 


v(3) = y. 

•'i.n y>,n- 


1 + 




GCKP84 is a general-purpose chemical kinetics 
program designed to solve a wide variety of problems 
(ref. 18). It uses the integration technique developed by 
Zeleznik and McBride (ref. 16) specifically to integrate 
chemical kinetic rate equations. Details of this 
integration technique are not yet available. 


DASCRU 

DASCRU is an explicit fourth-order Runge-Kutta 
process that differs from the standard Runge-Kutta 
method as follows: the standard Runge-Kutta method 
requires four derivative evaluations per step. The main 
disadvantage with its use is the difficulty of estimating 
the local truncation error at each step, since its formal 
expression is excessively complicated. The Runge-Kutta- 
Merson method requires an additional derivative 
computation that serves to determine the estimated local 
error. This technique uses the following equations to 
advance the solution from time /„_i to time t„ (ref. 19): 


v(l) = y. , + , 


In these equations 



A good estimate of the local error in computed the 
accepted value for is if the time step is 

sufficiently small (ref. 19). However, DASCRU uses the 
following test to ascertain convergence: 


ytl-y'il <0-5EPS for^f^) <10-3 


-yi.n -yi.n 
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(A34) 


<0.5EPS for >>(3^) >10-3 


where EPS is the local error tolerance. If any of the 
variables fails to satisfy the convergence criterion (eq. 
(A34)), the step length is halved and the step is retried. If 
all of the variables satisfy the convergence criterion (eq. 
(A34)), the step is accepted as successful and the solution 
is updated 


v(3) _ y . 
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Appendix B 

User-Supplied Parameters 


The codes (EPISODE, LSODE, CHEMEQ, CREKID, 
GCKP84, and DASCRU) examined in this study require 
the specification of several parameters in addition to the 
local error tolerance EPS and the elapsed time at which 
the integration is to be terminated. For each code, values 
for the user-supplied parameters that minimized the CPU 
time required by the code were obtained by a trial-and- 
error procedure. The parameters required by each code 
are discussed below. 

EPISODE 

EPISODE requires the user to specify the method flag 
MF, the error control to be performed lERROR, and the 
guess for the initial step length HO. For both test prob- 
lems all options (10, 11, 12, 13, 20, 21, 22, and 23) for 
MF were tried, and the stiff option with a user-supplied 
analytic evaluation of the complete Jacobian matrix 
(MF = 21) was found to be the most efficient. The type of 
error control to be performed is specified by the 
parameter lERROR, which has three possible values: 1, 
2, and 3. For lERROR = 1, the error control is absolute; 
and for lERROR = 2, it is relative. In the test problems 
examined in this study, the variables differ widely (see 
tables IV and V), so relative error control is appropriate. 
However, since some of the mole numbers had zero 
initial values, pure relative error control could not be 
used. The option lERROR = 3 was used instead. This is a 
semirelative error control. It is relative for variables that 
are initially nonzero. For a variable that is initially zero, 
the error control is absolute until the variable reaches 
unity in magnitude, when the control becomes relative. 
Since none of the mole numbers reaches a value of unity, 
the error control is always absolute for species with 
initially zero mole numbers. The optimum value for HO 
was found to be a function of the problem, the 
temperature method used, and the value for EPS. For 
test problem 2, values for HO close to the optimum value 
resulted in vastly increased CPU times. This behavior is 
illustrated in table VIII. Note that a decrease in HO from 
10-7 to 10-8 s has resulted in an order-of-magnitude 
increase in the CPU time. Although not shown here, the 
solution was also found to be adversely affected by a 
poor choice for HO. Also, some values for HO resulted in 
problems with solution instability. 

LSODE 

The user-specified parameters for LSODE include the 
method flag MF, the error control to be performed 
ITOL, and values for the local relative RTOL and 
absolute ATOL error tolerances. Both RTOL and ATOL 


can be specified either as (1) a scalar, so that the same 
error tolerance is used for all variables, or (2) an array, so 
that different values for the error tolerance can be used 
for different variables. The value of ITOL indicates 
whether RTOL and ATOL are scalars or arrays. ITOL 
has four possible values (1, 2, 3, and 4) which correspond 
to the types of RTOL and ATOL as follows: 

ITOL = 1 : scalar RTOL and scalar ATOL 
ITOL = 2: scalar RTOL and array ATOL 
ITOL = 3: array RTOL and scalar ATOL 
ITOL = 4: array RTOL and array ATOL 

The option ITOL = 2 (scalar RTOL and array ATOL) 
was used for the following reasons: since the same 
number of significant figures (given by RTOL) is 
acceptable for all components, RTOL was specified as a 
scalar. However, since the temperature has much larger 
values than the species mole numbers, ATOL can be 
much larger for the temperature than for the mole 
numbers. Hence ATOL was specified as an array. 

For both test problems, the options (10, 11, 12, 13, 20, 
21, 22, 23) for MF were tried and the method MF = 21 
(backward differentiation method and user-supplied 
analytic evaluation of the complete Jacobian matrix) was 
found to be the fastest. For species mole numbers, values 
for the absolute error tolerances that resulted in 
minimum CPU times were found to be a function of the 
problem, the temperature method used, and the value of 
the relative error tolerance EPS. For values of EPS given 
in tables XI and XII for test problems 1 and 2, 
respectively, a value of zero for the absolute error 
tolerance for the temperature (required by LSODE-B) 
resulted in minimum CPU times for both test problems. 

CHEMEQ 

CHEMEQ requires the user to specify the maximum 
number of corrector iterations ITMAX to be attempted 
before nonconvergence is declared and a smaller step 
length tried. The optimum value for ITMAX was found 
to be 5 for both test problems. Several attempts at 
increasing the efficiency of CHEMEQ were made. These 
included (1) replacement of the Newton-Raphson 
iteration procedure for the temperature (method A) with 
a single Newton-Raphson iteration and (2) limiting the 
temperature change per time step. If the temperature 
change during any time step exceeded the maximum 
permitted, the step was shortened accordingly and 
retried. 

Modification (1) resulted in less efficiency and more 
errors and was therefore abandoned. Test (2) was applied 
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(a) after each predictor step, (b) after each corrector step, 
and (c) after each predictor step and after each corrector 
step. For test problem 1 and temperature method A, 
option c was found to be the most efficient. 

Various values for the maximum temperature change 
permitted per time step were tried, and the optimum 
value was in the range 7 to 8 kelvins for test problem 1 . 
Although option c was found to be the best among the 
three options (a-c), it was only marginally faster than the 
basic algorithm (with no constraint on the temperature 
change per time step). Also, for test problem 2 with 
temperature method A, and for both test problems with 
temperature method B, none of the three options was 
more efficient than the basic algorithm. In addition, 
when the reaction rates were not continuously calculated 
but were updated only after an allowed temperature 
change AT, the basic method was found to be the most 
efficient. Hence modification 2 was also rejected. 

CREKID 

The user -supplied parameters to CREKID include the 
maximum number of corrector iterations ITMAX 
allowed before nonconvergence is declared; and the 
maximum temperature change DELT allowed before the 


thermodynamic properties /t,- and Cpj (/= 1,NS) and the 
rate constants kj l=AjT^j exp(-r/7)] and T_y, N_j, 
and k-j [=A^jT!^-j exp(-r_y/7)] are updated. Note 
that T_j, N-j, and Ar_y are calculated from the forward 
rate constants and the equilibrium constant. See Pratt 
and Wormeck (ref. 27) for details. The optimal value for 
ITMAX was in the range 5 to 7, depending on the test 
problem and the value of the relative error tolerance. For 
all runs presented in this report, a value of ITMAX = 10 
was used. The optimal value for DELT was in the range 1 
to 2 kelvins depending on the problem and on the value 
of the relative error tolerance. 

GCKP84 

Since details for this technique are not yet available, 
default values for all parameters were used. 

DASCRU 

This code requires the user to specify the guess for the 
initial step length HO. The optimum value for HO was 
found to a function of the problem, the temperature 
method used, and the error tolerance required of the 
solution. 
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Appendix C 

Approximation for Maximum Temperature Change Between Rate Constant Updates 


As discussed in the section Computational Tactics, 
significant increases in efficiency are realized by updating 
the forward rate constants kj [ = AjT^j exp ( - Tj/T)] and 
the backward rate constants 7’_y, N_y, and k^j 
[= A.ji^-J exp {-T-j/T)] only when the temperature 
change exceeds an amount AT that is problem dependent. 
In specifying a value for AT, care must be taken to avoid 
poor approximations in the resulting reaction rates 
because this leads to excessive computational work, as 
shown, for example, in table IX. To avoid the work 
associated with a trial-and-error search for an optimum 
value for AT, we now develop a simple approximation 
for it. 

Consider rate constant kj for reaction j. The exact 
expression for kj is given by 

kj = AjT^J exp (C 1 ) 



^/.approx 
^7, exact 


<EPS 


(C4) 


where the bars | | denote absolute value. 

The division of equation (C2) by equation (C3) and 
rearrangement give 



(C5) 


where tj=^T/T, the positive sign denotes increasing 
temperature, and the negative sign denotes decreasing 
temperature. 

There are four possible cases (other than the trivial one 
Tj=Nj = 0) as follows: 


where Aj is the preexponential constant, Tj is the 
activation temperature (= Ej/R, where Ej is the 
activation energy and R the universal gas constant), and 
T is the current value of the temperature. 

The idea behind using AT is that the work associated 
with computing kj and the reverse rate constants (T_y, 
A_y, and A’_y) from the forward rate constants and the 
equilibrium constant is eliminated for changes in T not 
greater than AT. Hence, if the temperature changes from 
its current value of Tby an amount not greater than AT, 
the rate constants kj and k-j are not reevaluated. The 
poorest approximation (A:y, approx) for the new rate 
constant is therefore given by equation (Cl) as follows: 

^7,approx~-^7^^^exp ^ 7 ^) (C2) 

The exact expression for kj is given by 

f^j,exact~ Aj(T±AT)^j exp (^ 3 ) 

where the plus sign denotes increasing temperature and 
the minus sign denotes decreasing temperature. 

We now estimate the maximum allowable AT by 
limiting the relative error in the resulting reaction rate to 
be no greater than EPS, the error tolerance required of 
the solution. 


(1) For increasing T 


(a) 1 


)>0 

(b) 1 


|<0 

(2) For decreasing T 

(a) 1 


1 >0 

(b) 1 


l<0 


Consider case a. For increasing T and positive 
(Tj/T +Nj), inequality (C4) becomes 


1 



For small EPS and ej-, exp(-EPS)sl -EPS, 
exp(€ 7 ’)sl + 67 ", and 1 /( 1 +€ 7 ) = 1 -ej-, so that equation 
(C6) becomes 
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exp exp (-£jN^)> exp (-EPS) 

or 

-er^^+N/^ ^-EPS 
or 

EPS 

for e^, EPS < < 1; 

^+Nj {Tj/T+n)^0 


Equation (C7) gives the maximum value of A 7 that can 
be used without resulting in a relative error greater than 
EPS in the forward reaction rate Rj for reaction j. For 
each forward and reverse reaction a value for the 
maximum allowable AT can be calculated by using 
equation (C7). Using the minimum of these values will 
ensure that no forward or reverse reaction rate will have a 
relative error greater than EPS. The minimum value for 
AT over all forward and reverse reaction rates is given by 


AT=min- 

J 


EPS T 


T 

J 


T^i 


or 


AT< 


EPS T 
T 

T ^ 


EPS T 


max 

J 




T 

J 



(C8) 


Similar analyses for the other three cases lead to similar 
expressions for AT. These four expressions can be 
replaced by a single expression given by 


AT 


EPS T 
T 

+A/. 


(C7) 


Equation (C8) provides a simple expression for the 
automatic evaluation of AT throughout the history of the 
problem. Every time the rate constants kj, T_j, and 
k-j are evaluated, which occurs only when the 
temperature change since the last update of the rate 
constants exceeds AT, a new value for AT can be 
calculated by using equation (C8). Using this equation 
therefore avoids the work associated with a trial-and- 
error search for the optimum value of AT. 
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